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EVALUATING  1  AND  2  DIMENSIONAL  MODELS  FOR  FLOODPLAIN 

INUNDATION 


CONTRACT  N68171-98-M-5830 


1.  Introduction 


1. 1  Objectives  of  the  study 

•  To  produce  a  state  of  the  art  numerical  model  of  a  large  reach  of  the  Missouri 
River,  capable  of  simulating  high  resolution  spatially/temporally  distributed 
results  of  water  depth  and  velocity  (TELEMAC  2D  modelling  system) 

•  To  investigate  model  behaviour  through  a  sensitivity  analysis. 

•  To  compare  model  predictions  to  midtiple  types  of  measured  data,  internal  to 
the  model  domain,  in  order  to  assess  model  performance  and  the  utility  of  the 
LANDS  AT  remote  sensing  validation  data. 

•  To  develop  a  model  application  using  LISFLOOD-FP  employing  SAR  remote 
sensing  data  on  a  reach  with  high  resolution  data  availability 

•  To  investigate  the  impact  of  varying  levels  of  bathymetric  data  and  mesh 
resolutions  on  the  LISFLOOD-FP  model  predictions. 

•  To  assess  the  current  R&D  needs  and  opportunities  in  eco-hydraulics 


1.2  Background  Information 

The  research  in  this  contract  focuses  on  achieving  improved  modelling  and  validation 
strategies  for  river  flood  inundation  and  the  within-channel  representation  of  flow 
processes.  This  style  of  research  and  investigation  is  of  considerable  importance  for 
applications  such  as  ensuring  appropriate  riverine  habitat  specifications,  eco- 
hydraulics  and  river  restoration  issues. 
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To  achieve  these  goals  the  research  programme  was  structured  in  the  following  terms: 

1 .  Background  and  key  features  of  the  study 

2.  Outline  of  the  TELEMAC  2D  modelling  system  for  river  hydraulic  modelling 

3.  The  Missouri  River  model:  Bathemetric  resolution  and  sensitivity  of  TELEMAC 

4.  The  Missouri  River  model:  Validating  TELEMAC  utilising  LANDSAT  TM  data 

5.  Outline  of  the  LISFLOOD  model  as  a  comparator  model  for  river  inundation  and 
within  channel  flow  process 

6.  Application  of  LISFLOOD  using  SAR  remotely  sensed  data  for  validation 
purposes 

7  Future  R&D  opportunities 


1.3  Key  Features  of  the  study 

•  The  application  of  a  state  of  the  art  two-dimensional  finite  element  code  for 
modelling  large  scale  fluvial  hydraulics  on  the  Missouri  River. 

•  The  use  of  a  high  resolution  model,  in  both  space  and  time,  along  with  wetting 
and  drying  algorithms  for  representing  moving  flow  field  boundaries  allows 
dynamic  inundation  predictions. 

•  For  the  first  time  model  predictions  are  validated  on  this  scale  in  both  time  and 
space  using  multiple  data  sources.  The  data  sources  are  internal  stage  data  at 
two  sites,  supplied  by  the  USAGE  Missouri  River  District  (MRD),  and  satellite 
imagery,  supplied  by  the  Remote  Sensing  and  Geographic  Information  Systems 
Center  at  USAGE  Cold  Regions  Research  and  Engineering  Laboratory 
(CRREL). 

•  The  strength  of  model  validation  is  graded  depending  on  the  data  source  and 
model  prediction. 

•  The  influence  of  bathymetric  data  on  the  model  predictions  is  studied  for  its 
potential  in  aiding  future  data  collection  strategy. 


2 


Anderson  and  Bates 


Flood  inundation  modelling 


2.  Outline  of  the  TELEMAC  2D  modelling  system  for  hydraulic  modelling 

The  TELEMAC  system  has  been  developed  by  the  Department  Laboratoire  National 
d'Hydraulique  (LNH)  at  Electricite  de  France,  Direction  des  Etudes  et  Recherches 
(EDF-DER).  LNH  was  formed  in  1946  to  undertake  studies  for  EDF's  hydroelectric 
power  projects  and  to  solve  the  hydraulic  engineering  problems  of  the  Maritime  Ports 
and  Waterways  Authority.  Since  1965  the  role  of  LNH  has  changed  somewhat.  About 
70%  of  LNH’s  work  is  now  directly  for  the  benefit  of  EDF,  split  between  hydraulics 
inside  power  station  machinery  and  the  study  of  environmental  problems  in  rivers, 
seas  and  the  atmosphere  brought  about  through  the  siting  of  such  plants.  The  rest  of 
LNH's  work  is  on  behalf  of  other  organisations. 

The  TELEMAC  system  therefore  follows  a  line  of  hydraulic  simulation  codes 
originating  from  LNH  for  the  study  of  environmental  problems,  such  as  flooding 
resulting  from  a  dam  break  or  thermal  emission  into  a  river  or  estuary.  TELEMAC 
was  developed  from  the  ULYSSE  code,  a  2D  finite  difference  system  (Ulysses  being 
the  father  of  Telemachus  in  Greek  mythology).  TELEMAC  is  now  in  its  third  four 
year  development  period,  each  period  having  a  budget  of  circa  US$16  million  (Hardy, 
1997).  TELEMAC  is  a  general  purpose  code  that  is  applicable  to  many  situations 
beyond  the  original  remit  of  its  design,  such  as  natural  floods  on  rivers. 

The  TELEMAC  system  is  a  series  of  computer  programs  utilising  finite  element 
techniques  for  simulating  hydraulic  situations.  The  system  includes  pre-  and  post¬ 
processing  components  and  offers  solutions  in  two-  and  three-dimensions.  It  also 
contains  facilities  for  sediment  and  contaminant  transport  (SUBIEF)  and  sand 
transport  (TSEF).  The  full  system  is  further  outlined  in  Figure  2.1.  All  the 
components  of  the  system  have  common  file  formats  and  are  written  in  the  high  level 
language  FORTRAN  90.  This  enables  easy  file  transfer  between  components  of  the 
system  and  allows  model  users  to  modify  parts  of  the  code  as  desired  for  specific 
applications. 

The  TELEMAC  system  contains  both  two-  and  three-dimensional  versions.  The 
three-dimensional  version  (TELEMAC-3D)  has  advantages  in  many  applieations 
where  vertical  velocity  variations  are  important  such  as  small  scale  and  oceanic 
studies  but  in  large  scale  fluvial  applications  depth  averaged  calculations  are 
adequate.  Indeed  the  equations  used  in  nearly  all  fluvial  hydraulic  models  derive  from 
the  depth  averaged  equations  known  as  the  Shallow  Water  or  St.  Venant  Equations, 
TELEMAC-2D  being  no  exception.  The  third  dimension,  though  initially  appealing, 
would  simply  add  to  the  computational  and  data  demands.  Descriptions  of 
TELEMAC-3D  can  be  found  in  Hervouet  et  al.  (1994),  Hervouet  and  Janin  (1994) 
and  Hervouet  and  Van  Haren  (1996).  TELEMAC-2D  is  however  used  throughout  the 
work  presented  here  and  is  now  discussed  in  more  detail. 
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2.1  TELEMAC-2D 

The  main  features  of  TELEMAC-2D  are  listed  by  Hervouet  et  al.  (1994)  as: 

•  structured  or  non  structured  meshing, 

•  use  of  Cartesian  or  spherical  co-ordinates, 

•  subcritical  and  supercritical  regimes  (with  hydraulic  jumps), 

•  various  momentum  source  terms  :  bottom  friction,  wind  stress  atmospheric 
pressure,  Coriolis  force, 

•  turbulence  modelling  (zero  equation  model,  k-epsilon  model), 

•  equation  on  temperature  or  a  substance  concentration, 

•  treatment  of  tidal  flats, 

•  many  types  of  boundary  conditions  including  free  slip  condition  and  incident 
waves. 

Many  of  these  are  of  minimal  importance  for  fluvial  applications  but  are  necessary  for 
the  diverse  nature  of  cases  that  TELEMAC-2D  can  be  applied  to.  A  wide  range  of  test 
cases  are  illustrated  in  the  TELEMAC-2D  (version  3.0)  validation  document  (Cooper, 
1996)  prepared  to  substantiate  the  explicit  claims  made  about  the  applicability  and 
accuracy  of  the  computer  code.  The  test  cases  illustrated  range  from  the  Western 
European  Coast  and  the  Mersey  estuary  UK,  to  flows  around  bridge  piers  and  over  a 
weir.  On  a  scale  more  relevant  to  this  study  dam  breaks,  flow  at  a  river  confluence 
and  the  simulation  of  a  flood  event  on  the  River  Culm,  UK  are  also  shown.  A  new 
version  of  TELEMAC-2D  is  released  annually.  All  of  the  simulation  in  this  report 
have  been  carried  out  using  version  3.0,  released  in  1995. 


2. 1. 1  Shallow  water  equations 

TELEMAC-2D  solves  the  shallow  water  equations  (SWE),  the  depth  averaged 
version  of  the  fully  three  dimensional  Navier  Stokes  equations  of  fluid  flow.  The 
SWE  require  that: 

•  the  flow  is  incompressible 

•  the  water  column  is  well  mixed,  so  that  there  are  no  significant  density  variations 
in  the  vertical 

•  vertical  accelerations  are  negligible  (hydrostatic  pressure  approximation) 

•  the  effective  turbulent  stresses  can  be  represented  by  an  eddy  viscosity  concept 
(the  Boussinesq  assumption) 

•  bed  stresses  can  be  modelled  using  a  linear  or  quadratic  friction  law. 

All  of  these  conditions  are  commonly  met  in  rivers,  estuaries  and  seas  making  the 
choice  of  these  equations  over  the  full  three  dimensional  Navier-Stokes  equations 
valid  for  the  applications  of  the  model  (Cooper,  1996). 

The  choice  of  formulation  of  the  SWE  used  in  TELEMAC-2D  is  not  obvious.  A 
conservative  form  would  seem  better  but  divisions  by  the  water  depth  are  needed  to 
produce  the  velocity  field,  hence  causing  problems  in  dry  areas  of  the  model  domain. 
Hence  a  non-conservative  form  is  preferred.  Moreover  numerical  stability  analysis 
also  favours  the  non-conservative  version  of  the  equations.  Using  finite  element 
methods  mass  conservation  can  be  ensured  with  non-conservative  equations.  Two 
versions  of  these  non-conservative  equations  have  been  developed,  the  celerity- 
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velocity  version  and  the  depth-velocity  version.  The  depth-velocity  version  are 
marginally  favoured  as  they  have  better  mass-conservation  properties  (Hervouet  and 
Janin,  1994).  These  equations  are  shown  below; 

— +  u.grad(h)  +  h.div(u)  =  0  2.1 

—  +  u.grad(u)  +  g  -  div(v.grad(u))  =  Sx  -  g  2.2 


^  A/  \  ^ 

— +  u.grad(u)  +  g  — - 


div(v.grad(v))  =  Sy  -  g 


dy 


2.3 


where; 

h  ;  water  depth 

u,  V  ;  velocity  components 

g  :  gravity 

Zf  ;  bottom  elevation. 

Sx,  Sy  ;Source/sink  terms  (bottom  friction,  wind,  etc.) 
V  ;  eddy  viscosity 


h,  u  &  V  are  the  unknown  variables. 

Equation  2.1  being  the  continuity  equation,  2.2  and  2.3  the  force-momentum 
equations.  Full  derivations  can  be  found  in  Norton  et  al  (1973). 


2.1.2  Solving  the  equations 

The  TELEMAC  system  uses  a  finite  element  methodology  to  solve  the  shallow  water 
equations  to  produce  values  of  water  depth  (h)  and  two  velocity  components  (u  and  v) 
at  all  points  in  the  model  domain.  To  achieve  this  the  domain  must  be  discretized  into 
a  grid  or  mesh,  usually  made  up  of  linear  triangles  (unstructured  grid)  for  flexibility. 
The  mesh  is  created  outside  TELEMAC  using  specialised  mesh  generation  software 
such  as  I-DEAS  or  BALMAT.  The  mesh  is  then  incorporated  with  the  topographic 
data  to  produce  a  geometry  file  that  is  then  used  in  the  TELEMAC-2D  simulation. 
The  mesh  consists  of  a  series  of  node  points,  the  vertices  of  the  triangles  with  a 
elevation  (z)  value,  where  the  equations  are  solved  and  the  triangles  themselves  are 
the  elements.  The  geometry  notation  used  is  shown  in  Figure  2.1. 

All  three  dependant  variables  (h,u  and  v)  are  defined  at  each  point  in  the  domain  in 
terms  of  linear  interpolation  functions  associated  with  the  value  attained  at  each  node. 
Linear  interpolation  functions  between  nodes  can  be  a  limitation  where  the  bed 
gradients  are  varying  rapidly  compared  to  the  distance  between  nodes.  Mesh 
refinement  can  however  be  used  to  overcome  this  problem  in  most  cases.  The  time 
discretization  is  semi-implicit  so  a  system  of  3N  simultaneous  algebraic  equations  is 
obtained,  where  N  is  the  number  of  nodes  in  the  domain  (Cooper,  1996).  Furthermore, 
TELEMAC-2D  employs  element  by  element  (EBE)  techniques  to  minimise  and 
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simplify  computation.  More  detailed  descriptions  of  finite  element  methodology  can 
be  found  in  Norton  et  al.  (1973)  and  Finder  and  Gray  (1977). 


k 

. . . rv: 

iy  ^ 

h:water  depth 

i;  :v: 

Z  :  free  surface 
elevation 

i 

Zf :  Bed  elevation 

Figure  2-1  The  geometry  notation  used  in  TELEMAC-2D  (after  Hervouet  and  Van 

Haren,  1995). 

Several  solution  algorithms  are  available  in  TELEMAC-2D  including  (Cooper,  1996): 

•  a  fractional  step  method  using  characteristics  for  the  advection  terms  and  a 
Galerkin  method  for  the  propagation  and  diffusion  terms  in  the  equations 

•  several  variants  of  the  Streamline  Upwind  Petrov  Galerkin  (SUPG)  method 
(Brooks  and  Hughes,  1982) 

•  a  hybrid  scheme  that  combines  the  characteristics  and  SUPG  approaches. 

Further  details  of  these  options  are  available  in  Hervouet  and  Van  Haren  (1995). 

In  the  simulations  presented  later  in  this  report  a  fractional  step  method  (Marchuk, 
1975)  is  used  where  the  advection  terms  are  solved  initially  followed  by  a  second  step 
where  the  propagation,  diffusion  and  source  terms  are  solved.  For  the  advection  step 
the  Method  of  Characteristics  is  used  for  the  momentum  equation  and  an  SUPG 
method  for  advection  of  h  in  the  continuity  equation,  ensuring  mass  conservation.  The 
use  of  the  semi-implicit  SUPG  method  for  the  continuity  equation  achieves 
unconditional  stability  (Hervouet  and  Janin,  1994).  The  second  step  of  the  calculation, 
propagation,  is  then  solved  using  a  conjugate  gradient  type  method. 

The  different  solution  techniques  do  however  produce  very  similar  results  when  used 
on  the  same  simulations  with  the  same  mesh.  For  example  Bates  et  al.  (1996)  show 
the  SUPG  and  hybrid  methods  used  on  a  flood  event  on  the  River  Culm,  UK,  showing 
very  similar  outflow  hydrographs  and  inundation  extent. 

The  Courant  number  is  used  as  a  measure  of  the  quality  of  the  numerical  solution.  It 
derives  from  explicit  numerical  methods  and  is  calculated  thus: 
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where  (Bates  et  al.  1996): 

C,  :  Courant  number 

t  :  time  step 

X  :  mesh  size 

u  :  flow  velocity. 

Explicit  methods  become  unstable  when  C^>\.  Although  TELEMAC-2D  is  implicit 
and  theoretically  not  subject  to  such  constraints  the  Courant  number  remains  a  good 
measure  of  the  quality  of  the  numerical  solution  (Hervouet  and  Janin,  1994;  Bates  et 
fl/.,1996).  Hervouet  and  Janin  (1994)  suggest  that  in  some  cases  simulations  can  be 
performed  with  values  of  the  Courant  number  of  up  to  50  but  this  is  not  advisable. 


2.1.3  Boundary  conditions 

There  are  two  main  types  of  boundary  condition  that  can  be  used,  solid  boundaries 
and  liquid  boundaries.  The  boundaries  described  here  are  round  the  sides  of  the  model 
domain  and  through  the  bed  of  the  model.  All  boundary  conditions  are  assigned  on  a 
node  by  node  basis.  Solid  boundaries  are  no  flux  (impermeable)  and  incorporate  a 
friction  factor  (Hervouet  and  Van  Haren,  1995).  Liquid  boundaries  allow  a  flux  across 
them.  They  are  more  difficult  to  deal  with  as  they  suppose  the  existence  of  a  fluid  area 
that  is  not  part  of  the  calculation  domain  but  can  however  influence  it.  This  influence 
is  described  through  the  boundary  condition.  There  are  four  types  of  liquid  boundary, 
entry  and  exit  with  supercritical  flow  (Froude  number  >1)  and  entry  and  exit  with 
subcritical  flow  (Froude  number<l).  Incident  waves  and  prescribed  flowrates  can  be 
incorporated  through  these  4  boundary  types.  Hervouet  and  Van  Haren  (1995) 
describe  these  boundary  conditions  more  fully. 


2. 1.4  Physical  parameter  options 

Physical  representation  of  parameters  is  essential  in  models  such  as  TELEMAC-2D  in 
order  to  apply  them  to  different  applications.  TELEMAC-2D  includes  numerous 
physical  parameters.  The  most  important  are  discussed  in  this  seetion. 


2. 1.4. 1  Bottom  friction 

There  are  six  options  in  TELEMAC-2D  for  representing  bed  friction,  these  are: 

•  No  friction 

•  Linear  friction 

•  Chezy’s  law 

•  Strickler’s  law 

•  Manning’s  law 

•  Nikuradse’s  law 
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The  applied  friction  coefficient  is  converted  to  force  terms  in  the  x  and  y  directions  at 
each  computational  node  which  are  then  fed  into  the  Shallow  Water  Equations  (in  the 
source  terms  Sx  and  Sy  -  see  section  2.1.1.,  equations  2.2  and  2.3).  In  reality  the  bed 
friction  force  is  a  quadratic  function  of  velocity  so  the  no  friction  and  linear  law  are 
rarely  used  in  practical  applications  (Hervouet  and  Van  Haren,  1996).  The  Chezy, 
Strickler  and  Manning  laws  are  all  closely  related  and  utilise  the  quadratic  function 
mentioned  and  are  all  described  more  fully  by  Hervouet  and  Van'  Haren  (1995). 
Nikuradse’s  law  calculates  a  Chezy  coefficient  from  the  water  depth  (h)  and  grain  size 
of  the  bed  material  which  is  then  converted  into  a  force  term.  Choice  of  law  is  not 
critical  as  they  are  very  closely  related  and  friction  coefficients  can  easily  be 
converted  between  them.  The  choice  rests  with  the  model  user. 

2.1. 4.2  Turbulence  representation 

Turbulence  is  one  of  the  major  unresolved  problems  in  physics.  Its  importance  in 
hydraulics  is  well  known  but  it’s  representation  in  hydraulic  models  is  a  different 
problem.  The  concern  in  hydraulic  modelling  is  how  turbulence  in  flow  affects  the 
mean  structure  of  the  flow.  Turbulence  modelling  is  a  scientific  discipline  in  itself  and 
the  reader  is  referred  to  Younis  (1996)  for  a  general  discussion  of  turbulence 
modelling  and  to  Hervouet  and  Van  Haren  (1996)  for  a  fuller  exposition  of  turbulence 
representation  in  TELEMAC-2D  than  is  possible  here.  In  TELEMAC-2D  version  3.0 
there  are  two  options  for  turbulence  representation: 

•  constant  viscosity 

•  k-epsilon  model. 

The  first  and  simplest  is  a  constant  viscosity  coefficient,  zero-equation  turbulence 
model.  The  term  VELOCITY  DIFFUSIVITY  is  used  to  set  the  molecular  viscosity, 
turbulent  viscosity  and  dispersion  throughout  the  model  domain.  In  turbulent  flows 
the  turbulent  viscosity  is  dominant  and  can  virtually  be  equated  to  the  velocity 
diffusivity  value.  A  constant  value  of  turbulent  viscosity  over  the  model  domain  is 
often  used  but  is  an  oversimplification. 

The  second  available  closure  is  the  k-epsilon  model  where  the  turbulent  viscosity  is 
expressed  as  a  function  of  the  turbulent  kinetic  energy  (k)  and  its  dissipation  rate  (e). 
This  two-equation  turbulence  model  is  overcomplicated  for  large  scale  applications 
and  is  computationally  intensive. 

A  third  turbulence  model,  the  single-equation  Elder’s  model,  has  just  been  introduced 
into  the  newest  version  of  TELEMAC-2D  (version  3.1).  This  model  separates 
longitudinal  and  transverse  dispersion  terms.  It  offers  improvements  on  the  previous 
two  methods  in  fluvial  applications  where  such  separation  is  important.  Unfortunately 
this  method  was  not  available  for  the  work  in  this  report. 


2. 1.4.3  Other  physical  parameters 

Several  other  physical  parameters  can  also  be  specified  in  TELEMAC-2D.  These 
include  wind  stress  on  the  water  surface,  atmospheric  pressure,  water  temperature  and 
water  density.  The  Coriolis  force  can  also  be  applied  when  modelling  large  areas. 
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2. 1.5  Wetting  and  drying  zones 

Applications  of  TELEMAC-2D  to  rivers  and  estuaries  generally  involve  wetting  and 
drying  of  areas  of  the  model  domain.  For  example  the  tidal  simulation  in  the  Mersey 
estuary  and  the  flood  event  on  the  River  Culm,  both  described  by  Cooper  (1996), 
involve  this  type  of  behaviour.  The  ability  of  the  model  to  deal  with  this  sort  of 
behaviour  is  therefore  of  vital  importance  but  is  problematic.  The  behaviour  in  the  dry 
zones,  where  divisions  by  the  water  depth  (h)  in  the  calculations,  can  cause  spurious 
terms  to  appear  as  h  tends  towards  zero. 

Two  solutions  to  this  problem  are  available  in  TELEMAC-2D,  both  deseribed  in  more 
detail  by  Hervouet  and  Van  Haren  (1995), 

•  solving  the  equations  everywhere  and  coping  with  the  spurious  terms. 

•  removing  the  dry  zones  from  the  computational  domain. 

The  first  is  the  simplest  but  corrections  must  be  applied  in  the  dry  zone.  In  dry  areas 
the  water  surface  gradient  becomes  that  of  the  bottom  topography  but  this  cannot  be 
allowed  to  act  as  a  driving  force  in  the  momentum  equation. 

The  second  method  removes  the  dry  zone  from  the  computational  domain  and  is  often 
called  the  "moving  boundary  technique".  In  TELEMAC-2D  this  is  achieved 
efficiently,  avoiding  the  need  to  redefine  the  mesh  at  each  time  step,  by  keeping  the 
elements  in  the  mesh  but  cancelling  their  existence  through  the  use  of  an  array  set  to  0 
for  dry  elements  and  1  for  all  others. 

Partially  dry  elements  are  another  important  area,  especially  where  the  elements  are 
large  such  as  in  the  Missouri  River  model.  These  are  coped  with  in  TELEMAC-2D  by 
a  sophisticated  method  that  allows  the  water  depth  to  go  to  zero  within  an  element. 
This  compares  favourably  with  the  usual  techniques  of  keeping  the  elements  fully  wet 
or  excluding  partially  dry  elements  (Figure  2.2). 


2.2  Post-processing  of  results 

After  the  computations  have  been  completed  the  numerical  results  must  be  converted 
to  a  more  user  friendly  visual  form.  This  is  done  using  the  graphics  visualisation 
software  RUBENS,  also  written  by  EDF-DER.  RUBENS  allows  quick,  easy  and 
flexible  representation  of  the  results  in  many  different  forms  such  as: 

•  mesh  plots 

•  vector  fields 

•  contour  lines 

•  coloured  surfaces 

•  space  profiles 

•  time  profiles 

RUBENS  also  allows  the  superposition  of  measurements  and  graphics  and  also  the 
manipulation  and  processing  of  the  visualised  data.  Most  of  the  plots  of  model  results 
in  this  report  are  created  using  the  RUBENS  software.  For  more  information  on 
RUBENS  the  reader  is  referred  to  Piro  (1993). 
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2.3  Summary  of  the  TELEMAC-2D  system 

TELEMAC-2D  is  a  high  resolution  space/time  distributed  hydraulic  model  that  solves 
the  Shallow  Water  Equations  for  fluid  flow  using  a  finite  element  methodology.  The 
model  can  be  used  in  a  wide  variety  of  scenarios  including  those  involving  wetting 
and  drying  fronts  within  the  model  domain.  Bed  friction  and  turbulence  are 
represented  in  the  model  through  the  use  of  physically  based  parameters. 

TELEMAC-2D  can  therefore  be  seen  to  be  well  specified  for  the  type  of  fluvial 
application  that  is  of  interest  in  this  study.  The  code  has  been  shown  to  work,  through 
the  report  of  Cooper  (1996),  in  a  wide  range  of  cases.  Success  in  any  individual 
situation  is  however  dependant  on  the  data  provided,  physical  and  numerical 
parameters  used.  How  these  important  factors  have  been  determined  in  the  Missouri 
River  model  case  is  the  subject  of  the  next  section 
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Figure  2-2  Partially  dry  elements:  (a)  the  real  situation,  (b)  element  remaining  fully  wet 
introducing  spurious  momentum  terms,  (c)  element  exclusion  method  and  (d) 
TELEMAC-2D  method  keeping  the  element  partially  dry  (after  Price,  1997). _ 
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3.  The  Missouri  River  model  :  Bathemetric  resolution  and  sensitivity  of 
TELEMAC 


This  section  looks  at  how  TELEMAC-2D  has  been  applied  to  the  Missouri  River 
between  Gavins  Point  Dam,  South  Dakota,  and  Maskell,  Nebraska.  The  theoretical 
aspects  of  the  model  were  examined  in  section  2  but  the  practical  details  of  applying 
the  model  to  this  specific  scenario  reach  are  vital  to  gain  an  understanding  of  the 
model  and  its  capabilities. 

The  reach  of  the  Missouri  River  being  used  for  this  modelling  study  is  that  from 
Gavins  Point  Dam,  South  Dakota,  to  the  gauging  station  at  Maskell,  Nebraska  (Figure 
3.1).  The  reach  covers  river  miles  811  to  776  making  it  35  miles  or  55  km  long.  The 
channel  varies  in  width  between  300m  and  1200m.  The  channel  slope  is  very  low, 
dropping  only  about  12  metres  along  its  55  km  length,  giving  a  gradient  of  0.02%. 
The  bed  material  in  this  channel  is  sand  which  is  fairly  mobile  but  the  channel  banks 
have  been  strengthened  or  stabilised  along  much  of  the  reach.  There  are  several 
islands  along  the  reach  several  of  which  are  permanent.  There  is  one  major  tributary, 
the  James  River,  that  joins  the  main  stem  at  river  mile  800  adjacent  to  the  James  River 
Island. 

The  flow  out  of  Gavins  Point  Dam  is  regulated  to  minimise  the  risk  of  flooding 
downstream,  hence  the  reach  being  modelled  is  very  unlikely  to  attain  an  out  of  bank 
conditions.  The  flow  in  the  James  River  is,  however,  naturally  variable.  The  model  of 
this  system  can  therefore  remain  as  channel  only.  River  flow  data  is  available  at 
several  points  along  the  reach  on  an  hourly  basis.  Table  3.1  shows  the  data  availability 
and  gauge  locations  (also  Figure  3.1). 


Table  3-1  Gauge  station  location  and  data  availability. 


Gauging  Station 

Distance  from  top  of  reach 
(km) 

Data  Available 

Gavins  Point  Dam,  SD. 

■■  6 

Flow  rate  ; 

Yankton,  SD. 

--  ....  7.5  . 

Gayville,  SD. 

■  ‘21.5-^ 

Maskell,  Neb. 

.  ,  ,  '  55 

Scotland,  SD. 

53  km  up,  James 
from  confluence 

Flow  rate  and  Stage 
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Figure  3-la)  Location  and  (b)  detail  of  the  study  reach  of  the  Missouri  River,  USA.  The 
modelled  reach  is  from  Gavins  Point  Dam  to  Maskell  and  includes  the  gauging  stations 
at  Yankton  and  Gayville. 
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The  flow  from  the  James  River  gauge  at  Scotland,  South  Dakota,  is  routed  using  a 
simple  one-dimensional  kinematic  wave  model  to  the  required  location  on  the  main 
model  input  boundary  when  necessary. 

The  hydraulic  model  used  in  this  study  is  TELEMAC-2D.  This  is  a  high  resolution 
space/time  distributed  hydraulic  model  using  the  finite  element  methodology.  The 
model  is  potentially  capable  of  fulfilling  all  the  objectives  of  this  study.  The  model  is 
described  in  detail  in  section  2. 


3. 1  Producing  the  finite  eiement  mesh 

The  finite  element  mesh  for  the  model  must  be  created  by  first  defining  the  boundary 
of  the  model  domain  and  secondly  discretizing  this  area  into  elements. 

The  boundary  of  the  2D  model  was  defined  by  digitising  round  the  edge  of  the  river 
as  represented  on  United  States  Department  of  the  Interior  7.5  minute  quadrangle 
series  maps.  The  use  of  wetting  and  drying  algorithms  in  the  model  enable  the  flow 
field  boundary  to  be  calculated  within  this  outer  boundary  thus  allowing  this  boundary 
to  be  fairly  loosely  defined.  The  digitised  boundary  was  then  converted  to  Universal 
Transverse  Mercator  (UTM)  metre  co-ordinates  for  compatibility  with  the  metric 
scale  used  by  TELEMAC-2D.  The  boundary  of  the  model  (with  TELEMAC-2D’s 
metric  scale),  which  includes  three  permanent  islands,  is  shown  in  Figure  3.2. 


A  mesh  of  triangles  or  quadrilaterals  must  be  made  within  the  model  domain  to 
facilitate  the  finite  element  solution  technique.  The  finite  element  meshes  used  with 
the  Missouri  model  were  generated  inside  the  prescribed  boundary  using  the  mesh 
generation  package  BALMAT.  The  meshes  created  were  made  of  near  equilateral 
triangles  in  order  to  increase  the  accuracy  and  minimise  mass  conservation  errors. 
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Two  meshes  have  been  created  of  this  reach.  The  two  are  shown  in  Figure  3.3  and 
their  attributes  compared  in  Table  3.2^ _ 

a) 

m 
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b) 
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12500 

0  10000  20000  30000  40000  m 

Figure  3-3  Finite  element  meshes  of  the  Missouri  River  from  Gavins  Point  Dam  to 
Maskell,  (a)  low  resolution  Mesh  1  and  (b)  high  resolution  Mesh  2. _ 

Mesh  1,  the  lower  resolution  mesh,  is  used  for  the  majority  of  the  simulations  in  this 
report.  It  should  be  assumed  this  has  been  used  when  looking  at  results  unless 
otherwise  stated.  Mesh  2  is  only  used  in  the  sensitivity  analysis  (section  4)  and  section 
6  where  its  performance  is  compared  to  that  of  mesh  1. 


Table  3-2  Comparison  of  the  attributes  of  the  two  meshes  used  to  model  the  Missouri 


River  in  this  report. 


Attribute 

Mesh  1 

Mesh  2 

Number  of  Nodes 

:  5969 

9213 

Number  of  Elements 

10437  7  :  ^ 

16567 

Average  Element  Size  (m^) 

100.77  A  ‘ 

'  79.98 
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3.2  Bathymetric  Data 

The  bathymetric  data  in  a  2D  hydraulic  model  is  one  of  the  most  important  factors  in 
attaining  high  quality  simulations.  The  quantity  and  quality  of  the  bathymetric  data  is 
therefore  of  utmost  significance. 

Along  this  reach  of  the  Missouri  River  the  bathymetry  was  obtained  by  the  MRD  in 
1995  by  echo  sounder  surveys  of  channel  cross  sections.  In  total  343  cross  sections 
were  available  along  the  55  km  reach  creating  an  average  spacing  of  just  over  150m. 
Each  cross  section  was  made  up  of  between  30  and  200  data  points,  depending  on 
section  length.  The  typical  spacing  between  data  points  along  a  section  was  6m.  All 
elevation  data  were  converted  to  metres.  Bathymetric  data  were  missing  on  the  cross 
sections  where  the  bed  was  above  the  water  surface,  such  as  on  sand  banks  and  mud 
flats.  In  these  regions  the  topographic  elevation  was  estimated  from  United  States 
Department  of  the  Interior  7.5  minute  quadrangle  maps,  except  on  the  permanent 
island  in  the  model  where  no  such  data  is  required. 

The  bathymetric  data  must  now  be  applied  to  the  finite  element  mesh.  This  is  the 
process  of  interpolating  between  topographic  data  points  and  assigning  elevation  (z) 
values  to  the  nodal  points  in  the  mesh.  This  produces  a  geometry  file  for  the  model 
simulation.  The  interpolation  is  usually  carried  out  using  STBTEL,  the  TELEMAC 
sub-program.  This  uses  a  quad-directional  interpolation  routine,  taking  the  nearest 
data  point  in  the  four  quadrants  around  each  node,  weighting  for  distance  from  the 
node  and  combining  to  give  the  nodal  elevation.  An  alternative  interpolation  routine 
has  been  written  at  Bristol  for  use  with  cross  sectional  data  of  river  channels.  This  is  a 
linear  interpolation  down  the  line  of  greatest  depth  between  sections  and  improves  the 
definition  of  this  line  in  such  cases  as  the  Missouri  River.  A  portion  of  the  models 
topography  is  shown  in  Figure  3.4. 
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Figure  3-4  An  example  portion  of  the  models  bathymetric  representation.  This  section  is 
immeadiately  upstream  of  Yankton. _ 
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3.3  Boundary  Condition  Specification 

The  boundary  conditions  of  the  model  are  very  important  to  the  simulations  run. 

There  are  several  boundaries  to  each  model  that  must  be  set  with  care  prior  to  the 
model  run.  These  boundaries  are  divided  into  two  types,  liquid  boundaries  and  solid 
boundaries. 

With  the  Missouri  model,  the  liquid  boundaries,  allowing  flow  across  them,  are  at  the 
top  and  bottom  ends  of  the  reach  and  on  the  James  River  tributary  inflow.  Model 
inflows  are  prescribed  as  flowrates,  at  Gavins  Point  Dam  and  from  the  James  River, 
and  outflows  as  water  surface  elevations,  at  Maskell.  This  produces  a  well  posed 
problem  for  fluvial  flows  according  to  the  theory  of  characteristics  (Hervouet  and  Van 
Haren,  1995).  All  measured  flowrates  were  converted  to  mVs  (cumecs)  and  stage 
values  to  metres. 

Solid  boundaries  allow  no  flux  across  them.  They  are  found  down  the  sides  of  the 
reach  and  at  the  bed.  The  side  boundaries  in  the  Missouri  model  are  set  as  slip 
boundaries,  allowing  a  velocity  along  them.  This  is  justified  by  the  size  of  the 
elements  making  flow  overly  restrained  in  areas  of  channel  constriction  when  the 
more  realistic  no-slip  boundaries  are  imposed.  The  flexible  boundary  of  the  flow  field 
removes  this  problem  for  large  portions  of  the  reach  where  such  an  argument  is 
irrelevant. 


3.4  Physical  Parameter  Specification 

The  two  most  important  physical  parameters  in  fluvial  hydraulic  models  are  generally 
agreed  to  be  bed  friction  and  turbulence  (Baird  and  Anderson,  1990;  Bates  et  al, 
1992).  The  theory  of  the  two  in  TELEMAC-2D  has  been  discussed  in  Section  2.  They 
are  the  only  two  considered  in  any  detail  here. 

Given  the  channel  only  nature  of  the  model  and  lack  of  additional  information  the  bed 
friction  and  turbulence  parameters  were  defined  as  constant  throughout  the  entire 
reach.  These  are  obviously  simplifications  but  should  be  adequate  as  a  first 
approximation.  The  bed  friction  is  always  defined  using  Manning’s  law  enabling  the 
standard  Manning’s  ’n’ measure  of  flow  resistance  to  be  used.  The  turbulence  was 
defined  using  the  zero-equation  velocity  diffusivity  representation. 


3.5  Sensitivity  Analysis 

The  first  stage  in  model  application  is  the  sensitivity  analysis.  Sensitivity  analysis  is  a 
widely  used  technique  in  hydrological  and  hydraulic  modelling  to  determine  the 
impact  of  changing  parameter  values  and/or  input  stresses  on  the  model  predictions. 
Virtually  every  modelling  study  involves  a  sensitivity  analysis  at  some  stage.  It  is  a 
technique  that  is  potentially  useful  in  model  formulation,  model  calibration  and  model 
verification  (McCuen,  1973). 
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3.3.1  Background  to  the  sensitivity  analysis 

Howes  and  Anderson  (1988)  suggest  that  a  sensitivity  analysis  can  be  used  to: 

•  Demonstrate  that  in  response  to  representative  variation  of  model  input  and 
parameter  values,  theoretically  realistic  model  behaviour  is  experienced. 

•  Illustrate  the  model  to  be  sufficiently  sensitive  to  represent  actual  variation  in  the 
prototype  system. 

•  Identify  those  model  parameters  or  inputs  to  which  the  model  is  most  sensitive. 

•  Assess  model  behaviour  without  recourse  to  comparisons  to  field  data. 

Perhaps  more  importantly  in  the  case  of  this  model  is  that  sensitivity  analysis  is  very 
useful  for  deciding  which  parameters  to  adjust  during  calibration.  Adjusting  the  most 
sensitive  parameters  will  produce  a  greater  improvements  in  the  predictions  for  small 
changes  in  parameter  values,  which  means  that  if  the  initial  values  were 
approximately  physically  realistic  then  the  calibrated  ones  should  be  as  well. 
Conversely  if  the  less  important  parameters  are  left  fixed  at  incorrect  values  then  the 
resulting  error  on  model  predictions  is  likely  to  be  small  (Troutman,  1985). 

Some  of  the  theoretical  considerations  of  sensitivity  analysis  are  now  discussed. 
Sensitivity  is  defined  as  the  rate  of  change  in  one  factor  with  respect  to  change  in 
another  factor.  Mathematically  it  can  be  derived  from  a  Taylor  series  expansion  of 
model  behaviour  and  after  discarding  the  non-linear  terms  the  linearized  sensitivity 
(S)  can  be  expressed  in  simple  terms  as  (McCuen,  1973) 

S  =  aF„/3F,  3.1 

where  Fq  the  model  output  prediction  and  R  is  the  model  input  parameter.  In  many 
instances  this  linear  approximation  of  sensitivity  is  adequate  to  describe  model 
behaviour  but  whether  it  is  appropriate  for  the  Missouri  model  is  unknown  at  present. 

There  are  two  methods  for  determining  the  value  of  the  model  sensitivity  which  are 
(Lane  et  ai,  1994). 

•  The  direct  differentiation  of  the  governing  equations.  This  has  been  demonstrated 
with  simple  models,  for  example  Beven  (1979)  uses  this  method  to  assess  the 
sensitivity  of  the  Penman-Monteith  evapotranspiration  equations  and  LaVenue  et 
al.  (1989)  use  it  to  help  estimate  travel  time  uncertainties  in  ground  water  flow. 
However  with  complex,  distributed  models  this  methodology  has  not  been 
developed  sufficiently  and  cannot  possibly  assess  the  entire  complex  response  of  a 
model  such  as  TELEMAC-2D. 

•  Factor  perturbation  approaches  can  evaluate  the  sensitivity  of  the  model  by 
incrementing  parameter  values  and  assessing  the  model  response.  This  is  the  most 
commonly  used  form  of  sensitivity  analysis. 

Given  the  problems  of  direct  differentiation  in  assessing  the  sensitivity  in  complex 
models  the  factor  perturbation  approach  is  the  only  one  considered  in  this  rest  of  this 
section. 

Assuming  the  overall  physically  based  validity  of  TELEMAC-2D  to  this  type  of 
application  the  sensitivity  analysis  of  the  Missouri  River  model  has  one  major 
function.  This  to  set  up  the  basis  for  the  following  calibration  procedures  by 
detenuining: 

•  tbe  relative  sensitivity  of  the  different  physical  parameters. 
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•  the  effect  of  mesh  resolution  on  this. 

Once  these  questions  have  been  resolved  the  application  and  calibration  of  the  model 
can  proceed. 

The  runs  carried  out  in  this  sensitivity  analysis  are  outlined  in  Table  3.3.  The  data 
from  these  simulations  is  enough  to  fulfil  the  two  aims  above.  As  discussed  in  the 
previous  section  only  two  physical  parameters  in  the  model  were  varied,  bed  friction 
and  turbulence  representation  (velocity  diffusivity),  both  of  which  were  evenly 
distributed  across  the  domain. 


Table  3-3  The  simulations  carried  out  for  the  sensitivity  analysis.  The  figures  in  bold  are 
the  parameter  values  held  constant  whilst  the  other  is  varied. 


Parameter  Varied 

Mesh  1 

Mesh  2 

Bed  Friction  (n) 

0.01,0.015,0.02, 0.025, 
0.03,0.035,0.04 

0.01,:  0.0 15,  0.02, 0.025, 
0.03,0.035, 0.04 

Turbulence  (Velocity 
Diffusivity  -  m^/s) 

'a,2,4-"- 

1,2,4 

The  range  of  bed  friction  was  chosen  to  represent  the  very  broad  range  for  this 
channel  type,  sand-bed  with  a  shallow  gradient  (0.02%  in  this  case)  (Table  3.4).  This 
full  range  is  used  here  as  no  data  are  available  to  make  a  further  judgement  at  this 
stage  and  uncertainties  in  the  overall  model  set-up  create  the  possibility  that  the 
values  needed  to  fit  the  results  to  the  observed  are  not  those  that  would  be  selected 
from  the  field.  The  velocity  diffusivity  values  are  harder  to  justify  physically  given 
their  vastly  oversimplified  application  to  the  flow.  The  values  chosen  are  similar  to 
those  used  in  other  applications  that  have  performed  well,  producing  realistic  output 
and  maintaining  model  stability. 


Table  3-4  Values  of  Manning’s  n  for  various  channel  types  (after  Bathurst  1988). 


Channel  Type 

Channel  Slope  (%) 

Manning  ‘n’  range 

Sand-bed 

<0.1 

0.01-0.04  ■ 

Gravel-bed 

0.05-0.5 

Boulder-bed 

IV  0.5-5  • 

'  '  0.03  -  0.2  '  "  - 

Steep  pool/fall 

>5 

'  0.1 -5  . 

The  parameter  space  was  sampled  using  simple  single  parameter  perturbation 
techniques.  The  flow  record  used  for  these  sensitivity  tests  is  that  from  around  the  6th 
June  1994  when  the  river  is  essentially  at  a  steady  state.  The  results  are  therefore  all 
steady  state  values  of  the  variables. 


3.5.2  Results  of  the  sensitivity  analysis 

The  results  of  the  parameter  changes  are  looked  at  on  the  inundated  area,  water  depth 
at  Gayville  and  Yankton,  velocity  at  Gayville  and  the  spatial  distribution  of 
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differences  of  water  depth  and  velocity  between  two  runs  over  the  whole  domain.  At 
this  stage  no  comparison  to  observed  data  is  made. 

Firstly  the  relative  sensitivity  of  the  two  parameters  on  both  meshes  is  looked  at. 
Figure  3.5  shows  the  impact  of  parameter  variation  on  the  total  percentage  inundation 
extent  in  the  domain.  It  can  clearly  be  seen  from  this  that  the  bed  friction  has  a  far 
greater  impact  than  the  velocity  diffusivity  over  the  ranges  tested. 


Figure  3.6  shows  the  impact  of  both  bed  friction  and  velocity  diffusivity  on  the  water 
surface  elevation  at  both  gauge  stations,  Yankton  and  Gayville,  for  both  mesh 
resolutions.  The  results  in  all  cases  show  that  the  bed  friction  has  a  far  greater 
influence  on  the  water  surface  elevation  than  the  velocity  diffusivity. 
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a)  Yankton 
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Figure  3-6  Sensitivity  analysis  plots  for  water  surface  against  bed  friction  at  (a)  Yankton 
and  (b)  Gayville  and  against  velocity  diffusivity  at  (c)  Yankton  and  (d)  Gayville. 

The  impact  of  parameter  variation  on  velocity  at  a  point,  in  this  case  at  Gayville,  is 
shown  in  Figure  3.7.  Bed  friction  is  again  dominant  but  the  response  is  not  linear. 
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Figure  3-7  Change  in  velocity  at  Gayville  as  (a)  bed  friction  and  (b)  velocity  diffusivity 
varies  for  both  high  and  low  resolution  meshes. 


From  the  above  three  diagrams  it  is  clear  that  the  bed  friction  is  clearly  the  most 
influential  parameter  having  a  far  greater  impact  than  velocity  diffusivity  on  all  the 
model  results.  The  mesh  resolution  has  a  consistent  impact  on  the  model  predictions. 
The  greater  the  resolution  of  the  mesh  the  lower  the  prediction  for  an  equivalent 
parameter  set  for  all  model  results.  This  is  a  slightly  unusual  result  and  shows  that  the 
model  structure  does  have  an  important  impact  on  the  model  predictions. 

Having  determined  that  the  bed  friction  is  by  far  the  dominant  parameter  in  the 
Missouri  River  model  further  analysis  was  carried  varying  only  on  this  parameter. 
Using  the  argument  of  Troutman  (1985),  the  velocity  diffusivity  can  only  introduce  a 
small  error  into  the  results  and  varying  it  produces  only  small  changes,  hence  it  can 
reasonably  be  overlooked.  The  difference  in  velocity  and  water  depth  were  calculated 
at  all  node  points  between  the  highest  friction  run  (n  =  0.04)  and  the  lowest  friction 
run  (n  =  0.01)  on  the  Mesh  1.  The  results  are  then  plotted  up  as  a  spatial  plot  of  the 
domain. 

The  water  depth  difference  plot  (Figure  3.8a)  clearly  shows  that  along  the  thalweg  the 
model  behaviour  is  expected  in  that  the  water  depth  is  greater  with  higher  friction 
(positive  difference  values).  On  areas  with  very  little  or  no  water,  such  as  sand  banks 
in  the  channel  or  out  of  channel  regions  in  the  model,  the  behaviour  is  very  much  less 
marked  with  only  very  small  changes  being  present,  as  would  be  expected.  There  are 
however  areas  in  the  model  domain  where  the  behaviour  looks  to  be  more  complex 
and  perhaps  unexpected. 
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The  velocity  difference  plot  (Figure  3,8b)  looks  to  more  complex  in  structure  than  the 
water  depth  one.  Along  the  thalweg  the  velocity  is  generally  lower  with  the  higher 
friction  (negative  difference  values)  but  in  the  shallower  regions  the  opposite  is  true. 
This  is  due  to  the  water  depth  increase  in  these  shallow  areas  facilitating  an  increase 
in  velocity  despite  the  increased  bed  friction.  There  are  however  many  areas  in  the 
model  where  the  behaviour  is  unexplained. 

Further  analysis  of  these  plots  has  been  carried  out  to  find  out  more  about  the 
complex  model  behaviour  at  all  the  nodes.  The  nodal  values  of  velocity  difference, 
water  depth  difference,  and  the  water  depth  for  the  low  friction  run  were  extracted 
from  TELEMAC-2D  results  files.  The  velocity  difference  was  then  plotted  against  the 
two  water  depth  variables  mentioned.  The  results  are  of  great  interest  (Figure  3.9)  and 
show  that  unexpected  model  behaviour  is  much  more  frequent  than  anticipated. 
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Figure  3.9a  shows  the  velocity  difference  plotted  against  the  water  depth  in  the  low 
friction  run  and  in  this  case,  froni  the  above  discussion,  it  would  be  expected  that  at 
all  nodes  but  those  of  very  shallow  depth  the  velocity  difference  should  be  negative. 
This  is  however  clearly  not  the  whole  story.  Although  the  majority  of  velocity 
differences  are  negative  there  are,  across  the  whole  depth  spectrum,  nodes  that  show 
the  opposite  behaviour.  At  very  low  (near  zero)  depths  the  majority  of  nodes  do  show 
an  increase  in  velocity  but  this  is  not  complete.  The  magnitude  of  the  velocity 
differences  seems  to  be  greater  at  smallish  depths  (around  Im)  possibly  caused  by  bed 
friction  having  a  greater  influence  on  flow  hydraulics  at  such  depths  than  at  greater 
depths.  Although  with  the  limited  number  of  points  at  greater  depths  this  could  just  be 
a  random  effect. 


a) 


Water  Depth  Difference  (m) 


Figure  3-9  Sensitivity  analysis  plots  showing  results  for  all  node  points  for  differences 
between  high  and  low  friction  runs  for  (a)  velocity  difference  against  water  depth  in  the 
low  friction  run  and  (b)  velocity  difference  against  water  depth  difference. _ 


Figure  3.9b  shows  the  velocity  difference  against  the  water  depth  difference.  The  first 
point  to  note  is  that  there  are  points  that  unexpectedly  show  a  decrease  in  water  depth 
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as  friction  increases.  All  these  points  also  show  a  decrease  in  velocity.  Between  water 
depth  changes  of  approximately  0  to  0.5m  there  is  predominantly  an  increase  in 
velocity  as  friction  increases,  probably  corresponding  to  the  nodes  with  shallow  water 
depths  initially.  At  greater  water  depth  differences  the  majority  of  the  behaviour  is  as 
expected,  decreasing  velocity  as  friction  increases  but  there  is  a  significant  proportion 
of  nodes  revealing  the  opposite,  though  not  to  the  same  magnitude. 

These  results  have  shown  the  complexity  of  the  model  response  to  very  simple 
changes  in  parameter  values.  Some  of  the  unusual  results  can  be  explained  by 
physical  reasoning,  others  are  perhaps  artefacts  generated  by  the  model  specification. 


3.5.3  Conclusions  of  the  sensitivity  analysis 

The  sensitivity  analysis  has  shown  that  the  bed  friction  is  by  far  the  most  influential 
parameter  on  all  outputs  in  this  model  application.  The  strategy  for  calibrating  the 
model  must  therefore  be  based  on  varying  this  parameter.  At  present  the  parameters 
are  distributed  evenly  over  the  entire  model  domain  but  although  this  is  perhaps 
simplistic  there  is  no  case,  or  data,  to  change  them  at  the  moment.  The  mesh 
resolution  has  a  consistent  influence  on  the  model  predictions  showing  the  model 
structure  is  an  important  factor  in  the  determining  the  model  predictions.  The 
sensitivity  analysis  also  highlighted  some  strange  and  unexpected  non-linear 
responses  in  some  areas  of  the  model  domain. 
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4.  Comparison  of  model  predictions  to  observed  data 
4. 1  Methods  of  Comparison 

The  data  available  for  this  study  is  stage  data  at  two  sites  within  the  model  domain 
which  can  be  considered  internal  but  1-D  and  numerical.  The  satellite  images  are  also 
internal  but  2-D  and  numerical.  The  satellite  imagery  is  therefore  the  more  powerful 
validation  data  type.  How  well  the  observed  and  predicted  data  match  shall  be  shown 
in  the  following  sections  and  at  the  end  this  idea  of  strength  of  validation  shall  be 
revisited  in  the  light  of  the  results.  Firstly  however,  the  methods  of  model  application 
used  in  this  study  and  elsewhere  are  reviewed. 

Usually  when  hydrological  models  are  applied  to  actual  scenarios  a  two  stage  process 
of  model  calibration  and  validation  is  carried  out.  Calibration  is  the  process  by  which 
parameter  values  are  varied  within  reasonable  ranges  until  the  differences  between 
observed  and  computed  values  are  minimised  (Konikow  and  Bredehoeft,  1992). 
Theoretically  speaking,  physically-based  models  should  not  need  calibrating.  The 
estimation  and  application  of  physically  realistic  parameter  values  for  the  model 
should  enable  the  model  to  perform  well  without  any  further  manipulation.  However, 
the  numerous  simplifications  involved  in  modelling  and  problems  in  parameter 
estimation  mean  that  such  a  notion  is  unrealistic.  Calibration  is  therefore  performed 
on  virtually  all  physically  based  hydrologic  models.  Following  the  calibration  phase 
the  model  must  then  be  validated.  In  most  cases  the  calibrated  parameter  values  are 
used  again  on  a  different  portion  of  the  flow  record  to  see  if  they  enable  a  good  match 
between  observed  and  predicted  variables.  If  they  do,  then  the  model  and  parameters 
are  considered  validated  and  can  be  used  with  confidence  for  prediction  of  future 
events.  If  not,  then  calibration  and  possibly  some  stage  of  model  re-formulation  must 
be  undertaken  until  the  model  can  be  validated.  It  should  be  noted  at  this  point  that 
validation  is  used  here  to  mean  that  a  reasonable  comparison  is  obtained  between 
observed  and  predicted  data  rather  than  that  the  model  has  been  shown  to  be  an 
accurate  representation  of  the  system. 

With  this  application  of  TELEMAC-2D,  the  highly  structured  methodology  for 
calibration  and  validation  outlined  above  is  not  entirely  appropriate  for  several 
reasons: 

•  this  study  is  for  research  of  model  behaviour, 

•  the  number  and  variety  of  data  sources  make  its  application  difficult, 

•  there  is  no  need  to  produce  a  model  capable  of  predicting  future  events  at  this 
stage, 

The  two  processes  are  therefore  not  used  as  distinct  entities  in  the  following 
investigation.  Instead  a  hybrid  of  the  two  forms  was  used  enabling  all  benefits  of  the 
calibration-validation  procedure  to  be  accrued  but  also  much  more.  The  model  was 
therefore  run  using  a  wide  range  of  bed  friction  values  over  several  sections  of  the 
flow  record.  This  enabled  the  calibrated  parameter  values  to  be  found  for  each 
comparison  data  set  (2  internal  stage  gauges  and  possibly  some  areas  on  a  satellite 
image).  Cross  comparison  of  parameter  values  and  the  assessment  of  errors  when  the 
model  was  calibrated  onto  one  data  source  allowed  an  assessment  of  the 
model/parameter  performance  more  complete  and  powerful  than  with  the  traditional 
calibration-validation  procedure. 
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The  time  period  being  simulated  for  this  comparison  is  around  6th  June  1994 
synchronous  to  the  LANDSAT-TM  satellite  image  of  the  area  that  will  be  used  in 
section  4.3.  The  flow  record  around  this  date  showed  very  little  (<2%)  variation  in  the 
data  set.  The  simulation  was  therefore  taken  as  a  simple  steady  state  computation, 
simplifying  the  analysis  of  results.  All  the  following  comparisons  are  therefore  done 
on  a  single  point  in  time  basis,  e.g.  a  single  water  level  for  each  point  in  the  domain 
for  each  parameter  set. 


4.2  Comparison  of  model  prediction  with  internal  stage  data 

Two  internal  gauges  for  water  level  are  available  for  the  reach  of  the  Missouri  River 
described  in  section  1.  The  gauges  are  at  Yankton,  7.5  km  down  the  reach,  and 
Gayville,  21.5  km  down  the  reach  (see  Table  3.1).  Hourly  stage  records  are  available 
for  these  gauges. 

The  sensitivity  of  the  model  prediction  of  water  depth  at  these  two  locations  has  been 
shown  to  be  dependant  predominantly  on  the  bed  friction  parameter.  This  is  the  only 
parameter  varied  in  the  following  analysis.  Other  parameters  would  have  only 
minimal  effect  but  increase  complexity  greatly.  The  sensitivity  of  water  depth  to 
friction  is  virtually  linear  at  both  internal  gauge  sites  (Figure  4-1)  and  in  the  expected 
direction,  i.e.  water  depth  increasing  as  friction  increases. 


Yankton 

Gayville 


Bed  Friction  (n) 

Figure  4-1  Linear  sensitivity  of  water  surface  elevation  to  bed  friction  at  the  two 
gauge  sites,  Yankton  and  Gayville.  _ 


This  finding  allows  the  use  of  McCuen’s  (1973)  linearized  sensitivity  equation  but  a 
simple  regression  equation  is  more  useful  for  applying  the  results  in  this  analysis  and 
is  therefore  used  here.  Equations  were  produced  for  the  model  behaviour  at  both  sites, 
Yankton  and  Gayville,  relating  the  error  in  water  level  prediction  to  the  bed  friction 
(Manning’s  ’n’)  parameter.  The  equations,  calculated  using  the  statistical  computing 
package  MINITAB,  can  be  written  thus; 

Yankton  Error  =  -0.410  +  34.72  n 
Gayville  Error  =  -0.976  +  31.93  n 
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Alternatively  they  can  relate  the  water  depth  to  the  bed  friction  thus: 

Yankton  Water  Depth  =  4.87  +  34.72  n 
Gayville  Water  Depth  =  2.91  +  31.93  n 

The  confidence  in  these  regression  lines  is  huge  with  tiny  residuals  (<0.01 5m)  in  all 
cases.  These  equations  can  therefore  be  used  as  a  very  powerful  tool  for  estimating 
model  predictions  for  ’n’  values  at  which  the  model  has  not  been  run,  although  it 
should  not  be  extrapolated  beyond  the  range  of  ’n’  values  used  to  make  the 
relationship  (i.e.  0.01  to  0.04).  At  present  this  is  only  shown  at  two  points  but  could 
perhaps  be  applied  at  every  node. 

It  can  easily  be  calculated  from  the  above  equations  that  the  ’n’ value  to  eliminate  the 
error  in  water  surface  elevation  at  Yankton  is  0.01 18  whereas  at  Gayville  it  is  0.0306. 
The  power  of  the  above  regression  equations  can  now  be  utilised  to  calculate  the  error 
in  the  other  observation  whilst  one  is  correct.  Whilst  the  stage  at  Yankton  is  predicted 
correctly  at  Gayville  the  prediction  is  0.60m  too  low.  Reversing  this,  whilst  Gayville 
is  predicted  correctly  the  prediction  at  Yankton  is  0.65m  too  high.  In  percentage  terms 
these  errors  are  18.26%  and  10.96%  of  the  expected  water  depth  at  the  two  points 
respectively. 

Taking  the  value  central  to  the  above  two  estimates  of  Manning’s  ’n’  should  enable  the 
approximate  calculation  of  the  joint  minimum  errors  in  the  predictions.  This  value  of 
’n’  is  0.0212  and  produces  an  error  at  Yankton  of  0.33m  and  at  Gayville  of  -0.30m. 
Slight  adjustment  of  ’n’  could  equalise  these  errors  but  would  be  irrelevant.  These 
errors  can  be  expressed  again  as  percentages  of  the  water  depth  such  that  the  error  at 
Yankton  is  5.89%  and  at  Gayville  is  8.37%.  Figure  4.2  shows  a  downstream  section 
of  the  comparison  between  observed  and  predicted  water  surface  values. 


Interestingly  the  value  of  rate  of  change  of  the  water  surface  (error  term  or  water 
depth)  is  different  at  the  two  gauging  stations.  Thus  a  fixed  change  in  the  bed  friction 
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produces  a  different  change  in  the  water  surface  elevation  at  the  two  stations. 
Investigation  into  how  this  type  of  behaviour  occurs  over  the  whole  reach  would  be 
useful  in  determining  more  about  model  behaviour  or  what  causes  this  phenomenon. 

The  comparison  between  observed  an  predicted  stage  values  has  shown  that  the  model 
is  capable  of  predicting  this  variable  to  about  0.3m  at  two  sites  simultaneously.  This  is 
a  good  result  given  the  simplicity  of  the  model  set  up.  Given  further  work  on  the 
parameterization  this  could  no  doubt  be  reduced  significantly. 


4.3  Comparison  of  model  inundation  predictions  with  satellite  imagery 

Many  types  of  satellite  imagery  have  been  used  in  studies  of  fluvial  and  coastal 
hydraulics.  The  wide  coverage,  ease  of  use  and  wealth  of  information  contained  in  the 
images  makes  them  ideal  for  many  uses  in  this  field.  Redfem  and  Williams  (1996) 
review  the  available  sensors  and  their  applications  in  (primarily  coastal)  hydraulics. 
Previously  a  lot  of  work  has  been  carried  out  on  the  remote  sensing  of  flood  extent 
because  of  the  immense  hazard  and  mitigation  costs  on  floodplain  developments. 
Rango  and  Salomonson  (1974)  show  that  the  areal  extent  of  flooding  can  be  mapped 
using  near-infrared  sensors  on  ERTS  1  (Earth  Resources  Technology  Satellite  -  later 
renamed  Landsat).  More  recently  Imhoff  et  al.  (1987)  use  SAR  (Synthetic  Aperture 
Radar)  and  Landsat  MSS  (Multi spectral  Scanner)  for  mapping  flood  extent  and 
damage  in  Bangladesh.  SAR  has  been  shown  to  very  accurate  for  flood  boundary 
delineation  by  Biggin  and  Blythe  (1996)  on  the  River  Thames  in  the  UK.  Satellite 
remote  sensing  of  the  1993  floods  on  the  Mississippi  and  Missouri  Rivers  in  the  USA 
with  both  SAR  (Brackenridge  et  al,  1994)  and  Landsat  TM  (Thematic  Mapper)  have 
further  shown  the  potential  of  such  methodologies. 

Satellite  data  has  also  been  used  to  measure  flood  stages  as  illustrated  by  Koblinsky  et 
al.  (1993)  using  the  U.S.  Navy’s  Geosat  radar  altimeter  on  the  Amazon  basin  and 
Brackenridge  et  al.  (1994)  using  SAR  data  and  topographic  maps  on  the  Mississippi. 
Both  methods  have  quite  large  errors  associated  with  them  at  present. 

Plumes  of  sediment  rich  water  can  also  be  identified  as  they  have  a  higher  reflectance 
than  clear  water  in  the  visible  region  of  the  spectrum  (Lillesand  and  Kiefer  1987). 
Brackenridge  et  al.  (1994)  use  this  factor  to  highlight  levee  breaches  along  the 
Mississippi  River  during  the  1993  floods,  the  breaches  acting  as  a  sediment  sources, 
and  explain  sediment  deposition  patterns.  Other  forms  of  pollution  and  thermal 
emissions  can  also  be  traced  (Lillesand  and  Kiefer,  1987). 

The  availability  of  satellite  imagery  for  the  modelled  reach  of  the  Missouri  River 
synchronous  to  the  flow  records  has  enabled  several  fundamental  research  objectives 
to  be  addressed.  Data  from  satellite  images  is  ideal  for  use  with  2-D  fluvial  hydraulic 
models  being  of  a  scale  commensurate  to  the  model  resolution  and  being  widely 
spatially  distributed.  In  this  study  such  data  has  the  potential  to  be  used  for: 

•  validation  of  inundation  extent  predictions, 

•  identification  of  regions  of  error, 

•  as  calibration  data, 

•  and  for  process  validation. 
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Remotely  sensed  data  has  the  potential  to  be  used  more  widely  in  hydraulic  modelling 
for  topographic  and  physical  parameterization  (Bates  et  al,  1997)  but  these  are 
ongoing  research  themes. 

Two  sites  at  which  to  analyse  model  inundation  predictions  were  chosen.  These  were 
selected  specifically  because  of  their  complex  bathymetry  and  topography,  involving 
mud  flats,  permanent  and  temporary  islands.  The  location  of  the  two  sites  are  shown 
on  Figure  4-3.  Side  by  side  comparisons  of  the  observed  and  predicted  inundation  are 
to  be  made  using  different  friction  values  at  both  sites  followed  by  an  overlay 
showing  the  geo-referencing  that  can  be  done  to  the  results. 


Figure  4-3:  Locations  of  the  two  sites  to  be  used  in  the  comparison  of  model  predictions 
with  satellite  imagery. 

The  side  by  side  comparisons  of  the  model  predictions  against  the  Landsat  TM  image, 
using  bed  friction  values  of  0.01,  0.02,  0.03  and  0.04,  are  shown  in  Figure  4-4  and 
Figure  4-5  for  the  two  sites.  These  figures  plot  the  satellite  image  against  the 
predicted  flow  field  boundary  as  this  is  the  model  result  that  can  be  directly  compared 
to  the  image.  From  these  it  can  be  seen  that  the  model  predictions  are  generally  good, 
producing  a  close  spatial  match  against  the  observed  data.  This  allows  calibration  of 
the  model  to  be  carried  out  with  such  data  by  varying  bed  friction  until  there  is  a  close 
spatial  match  between  observed  and  predicted.  This  variable  behaviour  is  caused  by 
the  gradient  of  the  topography  around  the  water  surface.  If  the  topography  has  a 
shallow  gradient  then  changes  in  the  bed  friction,  changing  the  water  depth,  causes 
changes  in  the  inundated  areas.  Where  the  topographic  gradient  is  steeper  then  the 
same  changes  in  water  depth  shall  have  virtually  no  impact  on  inundated  area.  A 
problem  arises  in  this  specific  case,  as  far  as  calibrating  the  model  on  this  data, 
because  the  topography  above  the  water  surface  was  only  estimated  from  US 
Department  of  the  Interior  maps  and  is  therefore  probably  not  of  the  accuracy 
required  to  enable  much  confidence  to  be  placed  in  the  predictions  of  inundated  area 
on  islands. 
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Figure  4-4:  Comparison  of  Landsat  TM  image  and  4  alternative  model  predictions  with 
different  bed  friction  (Manning’s  ‘n’)  values  at  site  1. 


Figure  4-5:  Comparison  of  Landsat  TM  image  and  4  alternative  model  predictions  with 
different  bed  friction  (Manning’s  ‘n’)  values  at  site  2. 

By  geo-referencing  the  model  results  and  satellite  image  an  overlay  of  the  two  can  be 
made.  An  example  of  such  is  shown  in  Figure  4-6  for  a  section  of  this  reach  including 
site  2.  This  detailed  comparison  between  observed  and  predicted  flow  field  boundary 
enables  problem  areas  to  be  pinpointed  and  further  model  development  to  proceed  by 
further  data  collection  or  adjustments  to  the  model  in  these  areas. 


Figure  4-6  An  overlay  of  the  model  prediction  of  flow  field  boundary  onto  the 
Landsat  TM  image  for  a  region  around  Goat  Island  (including  site  2). _ 


4.4  Summary 

These  first  simulated  flow  conditions  have  shown  both  the  excellent  predictions  of  the 
2-D  hydraulic  model  and  the  utility  of  the  data  sources  to  validate  its  performance. 
The  internal  stage  data  was  matched  at  single  points  but  not  so  closely  at  both  points 
simultaneously.  The  inundation  extent  was  validated  against  the  satellite  images  very 
well.  Some  areas  showed  up  problems  in  the  models  bathymetry/topography.  Given 
possible  future  refinements  in  the  model  set-up  the  internal  stage  data  should  be  able 
to  be  matched  simultaneously  at  both  sites  quite  easily.  Matching  all  the  flow  field 
boundaries  may  be  more  difficult  given  their  dependence  on  the  model  bathymetry. 
This  is  however  expected  given  the  relative  strength  of  the  two  validation  types  as 
discussed  earlier. 
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5.  Outline  of  the  LISFLOOD  model  as  a  comparator  model  for  river 
inundation  and  within  channel  flow  process 

To  extend  the  above  conclusions  the  TELEMAC-2D  model  was  compared  to  an 
alternative  simple  inundation  prediction  code,  LISFLOOD-FP,  for  a  UK  site  where 
high  resolution  inundation  extent  data  was  available  from  SAR  imagery.  The 
LISLOOD-FP  model  is  described  in  this  section,  whilst  the  high  resolution  data  and 
model  inter-comparison  are  described  in  Section  6.  Finally,  the  impact  of  varying 
levels  of  bathymetric  data  provision  and  mesh  resolution  on  the  LISFLOOD-FP 
model  is  discussed  in  Section  7. 


This  LISFLOOD-FP  model  is  described  fully  in  Bates  and  De  Roo  (2000),  but  the 
salient  points  are  reproduced  here  along  with  improvements  made  to  the  model 
subsequent  to  the  original  paper.  Channel  flow  is  handled  using  a  one-dimensional 
approach  that  is  capable  of  capturing  the  downstream  propagation  of  a  floodwave  and 
the  response  of  flow  to  free  surface  slope,  which  can  be  described  in  terms  of 
continuity  and  momentum  equations  as: 
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Q  is  the  volumetric  flow  rate  in  the  channel,  A  the  cross  sectional  area  of  the  flow,  q 
the  flow  into  the  channel  from  other  sources  (i.e.  from  the  floodplain  or  possibly 
tributary  channels).  So  the  down-slope  of  the  bed,  n  Manning’s  coefficient  of  friction, 
P  the  wetted  perimeter  of  the  flow,  and  h  the  flow  depth.  In  this  case,  the  channel  is 
assumed  to  be  wide  and  shallow,  so  the  wetted  perimeter  is  approximated  by  the 
channel  width.  The  term  in  brackets  is  the  diffusion  term,  which  forces  the  flow  to 
respond  to  both  the  bed  slope  and  the  free  surface  slope,  and  can  be  switched  on  and 
off  in  the  model,  to  enable  both  kinematic  and  diffusive  wave  approximations  to  be 
tested.  With  the  diffusion  term  switched  off,  equation  (5.2)  can  be  solved  for  the  flow 
Q  in  terms  of  the  cross  section  A,  and  hence  a  partial  differential  equation  in  A  is 
derived  from  equation  (5.1).  Usually,  Q  is  chosen  as  the  dependent  variable  (Chow, 
1988,  p  296),  as  it  results  in  smaller  relative  errors  in  the  estimation  of  discharge.  For 
this  model,  however,  we  are  primarily  interested  in  water  levels  (which  dictate  the 
flood  extent),  so  the  cross  sectional  area  A  is  used  as  the  dependent  variable.  It  also 
simplifies  the  inclusion  of  the  diffusion  term.  An  explicit  non-linear  finite  difference 
system  in  A  is  then  solved  using  the  Newton-Raphson  technique,  rather  than  the 
linearised  scheme  used  in  the  original  model.  The  diffusion  term  can  be  included  in  an 
explicit  fashion  simply  by  modifying  the  bed  slope.  So,  to  include  the  depth  gradient 
term.  This  approach  can  cause  instability  and  the  development  of  saw-tooth 
oscillations  in  the  solution,  but  these  are  easily  countered  by  the  use  of  a  smaller  time 
step.  If  this  is  the  case,  the  time  steps  for  solving  for  channel  and  floodplain  flows  can 
be  effectively  de-coupled  by  using  a  number  of  sub-iterations  for  the  channel  flow. 

A  flow  rate  is  imposed  at  the  upstream  end  of  the  reach,  which  for  the  kinematic  wave 
model  is  sufficient  as  a  boundary  condition,  as  wave  effects  can  only  propagate 
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downstream,  and  no  backwater  effects  need  to  be  taken  into  account.  If  the  diffusion 
term  is  included,  some  downstream  boundary  condition  is  required  to  close  the 
solution,  as  backwater  effects  are  taken  into  account.  This  can  either  be  a  stage 
reading,  or  a  zero  free  surface  slope  condition  can  be  imposed,  which  leaves  the  depth 
at  the  downstream  boundary  free  to  vary,  but  prevents  the  solution  developing  a  draw¬ 
down  or  draw-up  curve.  (This  option  will  be  referred  to  as  a  free  downstream 
boundary  condition.) 

The  channel  parameters  required  to  run  the  model  are  its  width,  bed  slope,  depth  (for 
linking  to  floodplain  flows)  and  Manning’s  n  value.  Width  and  depth  are  assumed  to 
be  uniform  along  the  reach,  their  values  assuming  the  average  values  taken  from 
channel  surveys.  A  uniform  bed  slope  is  calculated  from  the  DEM  (which  is  assumed 
to  be  too  coarse  to  contain  any  explicit  channel  information),  by  linear  regression 
along  the  line  of  the  channel,  which  is  defined  by  a  series  of  vectors  derived  from 
large  scale  maps  of  the  reach.  It  will  be  possible  to  derive  such  channel  parameters 
from  high  resolution  DEMs  automatically.  The  remaining  friction  coefficient  is  left  as 
a  calibration  parameter. 

Floodplain  flows  are  similarly  described  in  terms  of  continuity  and  momentum 
equations,  discretized  over  a  grid  of  square  cells,  and  two  options  exist  in  the  model 
for  the  treatment  of  floodplain  flows.  Most  simply,  we  can  assume  that  the  flow 
between  two  cells  is  simply  a  function  of  the  free  surface  height  difference  between 
those  cells  (Estrela  and  Quintas,  1994): 
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where  h‘’^  is  the  water  free  surface  height  at  the  node  (i,j).  Ax  and  Ay  are  the  cell 
dimensions,  n  is  the  Manning’s  friction  coefficient  for  the  floodplain,  and  Qx  and  Qy 
describe  the  volumetric  flow  rates  between  floodplain  cells.  Qy  is  defined  analagously 
to  equation  (5.4).  The  flow  depth,  hnow,  represents  the  depth  through  which  water  can 
flow  between  two  cells,  and  is  defined  as  the  difference  between  the  highest  water 
free  surface  in  the  two  cells  and  the  highest  bed  elevation  (this  definition  has  been 
found  to  give  sensible  results  for  both  wetting  cells  and  for  flows  linking  floodplain 
and  channel  cells.)  The  second  option  is  to  discretise  the  diffusive  wave  equation  over 
the  grid: 


V.J 


i.iX 


hlUh'-'-'-h*-' 


Ax 


Ay 


^  h  i-‘-j  _  U  i'-i  f  Vi  '-H  _  h  ‘■j+' 


sl/4 


Ax 


-H 


2Ay 


(5.5) 


The  two  approaches  are  subtly  different:  in  the  diffusive  approach,  the  x,y 
components  of  the  flow  are  linked,  whereas  in  the  cellular  approach,  the  flow  between 
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cells  is  solely  a  function  of  the  component  of  free  surface  gradient  in  that  direction.  A 
possible  criticism  of  the  cellular  approach  is  that  it  fails  to  reproduce  some  (intuitively 
correct)  features  of  floodplain  flows.  For  example,  flow  may  not  be  parallel  to  the  free 
surface  gradient,  depending  on  the  orientation  of  the  free  surface  slope  and  the  model 
grid,  but  the  two  approaches  agree  when  the  slope  is  parallel  to  one  of  the  grid  axes. 
The  differences  between  the  two  models  can  be  derived  analytically  for  a  free  surface 
slope  of  unit  magnitude  as  a  function  of  direction.  The  flow  vectors  differ  in 
magnitude  by  a  mean  value  of  -20%  and  in  angle  by  -10°.  While  these  deviations 
may  be  partially  compensated  for  in  the  friction  calibration  process  (the  cellular 
approximation  predicts  larger  flows  than  the  diffusive  approximation),  the  effect  on 
the  bulk  flow  behaviour  of  the  model  is  unclear,  and  so  both  approximations  are 
tested  in  this  study.  Whichever  approximation  is  adopted,  an  explicit  scheme  is  used: 
floodplain  flows  are  calculated  first  using  equation  (5.4  or  5.5),  then  the  water  depths 
on  the  floodplain  are  updated  using  (5.3). 

Equations  (5.4  or  5.5)  are  also  used  to  calculate  flows  between  floodplain  and  channel 
cells,  allowing  floodplain  cells  depths  to  be  updated  using  (5.3)  in  response  to  flow 
from  the  channel.  These  flows  are  also  used  as  the  source  term  in  (5.1),  effecting  the 
linkage  of  channel  and  floodplain  flows.  Thus  only  mass  transfer  between  channel 
and  floodplain  is  represented  in  the  model ,  and  this  is  assumed  to  be  dependent  only 
on  relative  water  surface  elevations.  While  this  neglects  effects  such  as  channel- 
floodplain  momentum  transfer  and  the  effects  of  advection  and  secondary  circulation 
on  mass  transfer,  it  is  the  simplest  approach  to  the  coupling  problem  and  should 
reproduce  some  of  the  behaviour  of  the  real  system. 

Comparison  of  results  from  this  simple  code  to  those  obtained  from  the  TELEMAC- 
2D  model  have  the  potential  to  provide  insights  into  both  the  physics  of  flood 
inundation  and  the  nature  of  the  modelling  process.  In  particular,  they  may  allow  one 
to  quantify  the  level  of  data  provision  necessary  to  achieve  a  particular  prediction 
accuracy  and  allow  one  to  more  fully  assess  the  ability  of  satellite  derived  inundation 
data  to  discriminate  between  competing  model  formulations. 


6.  Application  of  LISFLOOD  using  SAR  data  for  validation 

6. 1  Test  site  and  validation  data 

A  4  km  long  reach  of  the  upper  river  Thames,  UK,  has  been  selected  for  comparison 
of  the  two  models.  For  this  site  a  50m  resolution,  25cm  precision  airphoto  DEM  is 
available,  along  with  ground-surveyed  channel  cross  sections.  The  DEM  has  been 
modified  by  the  inclusion  of  a  dyke  (identified  in  maps  of  the  reach)  which  runs  for 
500m  along  the  north  side  of  the  channel  at  the  upstream  end  of  the  reach.  This  was 
found  necessary  to  constrain  the  flow  in  this  region,  and  if  omitted  causes  a  gross 
overestimation  of  flood  extent  in  the  upstream  part  of  the  model.  A  gauging  station  is 
located  at  the  upstream  end  of  the  reach,  which  can  provide  a  boundary  condition  for 
the  model.  The  floodplain  environment  is  entirely  agricultural,  being  mainly  made  up 
of  meadow  and  rough  pasture,  which  should  make  for  an  easier  calibration  j^roblem 
with  respect  to  floodplain  friction.  Bankful  discharge  is  estimated  at  40m^s'  .  A  finite 
element  mesh  for  the  TELEMAC-2D  model  has  been  constructed  using  a  streamline 
curvature  dependent  mesh  generator  (Horritt,  2000b)  that  ensures  depth  prediction 
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errors  are  independent  of  channel  sinuosity.  Element  sizes  vary  from  <50m  near  the 
channel  to  -lOOrn  in  the  shoreline  region.  The  mesh,  along  with  the  topography  and 
channel  for  the  LISFLOOD-FP  model,  is  shown  in  figure  6.1.  Steady  flow  is  assumed 
throughout  this  study  (peak  inflow  of  73  m^s’’),  as  given  the  short  length  of  the  reach 
and  a  typical  kinematic  wave  propagation  velocity  (1.5  ms'^),  we  would  expect  the 
reach  to  respond  to  changes  in  inflow  in  <1  hour.  Even  taking  into  account 
propagation  times  for  the  wetting  front  (measured  at  ~3  hours),  this  is  still  far  quicker 
than  significant  changes  in  inflow  (the  hydrograph  peak  is  ~40  hours  in  duration). 

Validation  data  are  also  available  in  the  form  of  satellite  imagery  (from  the  ERS-1 
SAR  sensor)  of  a  l-in-5  year  flood  event  on  the  reach.  This  has  been  processed  using 
a  statistical  active  contour  technique  (Horritt,  1999,  Horritt  et  al,  2000)  to  extract  the 
flood  shoreline,  which  can  then  also  be  used  to  form  a  raster  map  of  the  inundation 
state.  This  technique  has  been  found  to  delineate  the  shoreline  with  a  mean  location 
error  of  ~50m  (Horritt  et  al,  2000)  when  compared  with  airphoto  data,  which  is 
probably  adequate  for  this  study,  as  it  commensurate  with  the  DEM  resolution. 

6.2  Model  testing 

Friction  coefficients  for  the  channel  and  floodplain  remain  unconstrained  for  this 
problem,  and  are  therefore  treated  as  calibration  coefficients  (Bates  et  al,  1998, 
Horritt,  2000a).  While  this  is  likely  to  be  a  major  source  of  model  error,  and  will 
certainly  cloud  the  issue  of  model  comparison,  the  lack  of  alternative  techniques  for 
parameterising  friction  means  that  calibration  is  currently  the  only  way  forward.  It 
also  offers  the  opportunity  of  exploring  the  effects  of  friction  parameterisation  on  the 
modelling  strategies. 

Firstly,  we  define  the  extent  and  dimensionality  of  the  parameter  space  for  the 
calibration  problem.  We  assume  only  two  friction  classes,  one  for  the  channel  and  one 
for  the  floodplain,  giving  a  two  dimensional  problem.  Manning’s  n  values  for  the 
channel  range  from  0.01  to  0.05,  equivalent  to  values  quoted  for  concrete  lined 
straight  channels  and  winding  natural  channels  with  vegetation  and  pools, 
respectively  (Chow,  1988,  p35).  Values  for  the  floodplain  range  from  0.02  to  0.10  or 
0.12  (depending  on  the  model  used,  see  discussion  below),  equivalent  to  a  surface 
somewhat  smoother  than  pasture  (n=0.035)  to  dense  trees.  This  enables  the  full  range 
of  possible  frictional  values  to  be  explored.  For  this  reach  (winding  channel 
surrounded  mainly  by  pasture  with  hedgerows),  we  would  expect  the  optimum 
calibration  to  occur  approximately  in  the  centre  of  the  parameter  space,  but  as  the 
friction  calibration  is  also  partly  used  to  compensate  for  poorly  represented  processes 
in  the  model,  the  optimum  may  be  shifted.  Exploring  the  entire  parameter  space  is  a 
computationally  intensive  process,  so  a  more  pragmatic  approach  is  adopted,  instead 
aiming  to  explore  only  sparsely  the  full  space,  but  focussing  more  attention  on  the 
area  in  the  region  of  the  optimum  calibration. 

Before  calibration  can  be  performed,  we  first  also  need  to  define  some  measure  of  fit 
between  the  observed  and  predicted  flood  extent,  as  it  is  this  measure  that  will  be 
optimised  by  the  calibration  process.  We  use  an  area  based  measure,  the  area  correctly 
predicted  as  either  wet  or  dry  by  the  model,  which  is  corrected  for  bias  which  may  be 
introduced  by  the  area  occupied  by  the  flood.  For  example,  for  a  small  flood  in  a 
large,  mostly  dry  domain,  even  a  relatively  poor  prediction  of  flood  extent  may  give 
apparently  good  results  in  terms  of  the  area  correctly  predicted  by  the  model  (i.e.  the 
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Figure  6.1:  Finite  element  mesh  and  airphoto  topography  for  the  Thames  text 

site. 
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large  dry  area).  This  can  be  rectified  by  ignoring  an  appropriate  portion  of  the  dry 
area,  to  make  it  equal  in  size  to  the  flood.  In  this  case,  only  15.1%  of  the  domain  is 
flooded  and  84.9%  is  dry,  and  so  in  the  calculation  of  correct  area,  69.8%  (84.9- 
15.1%)  is  disregarded  when  the  model  predictions  and  satellite  data  are  compared,  in 
order  to  give  approximately  equal  flooded/unflooded  areas.  Comparing  results  from 
the  two  models  and  the  satellite  data  presents  a  little  difficulty,  due  to  the  different 
scales  and  representations  of  the  depth  field  used  in  the  models.  In  this  case,  the 
results  from  the  raster  model  are  sampled  onto  a  12.5m  grid  (the  resolution  of  the 
satellite  data),  and  then  can  be  compared  directly  on  a  pixel  by  pixel  basis.  The  finite 
element  results  are  also  sampled  onto  a  12.5m  grid,  the  depth  field  being  linearly 
interpolated  across  each  element.  The  treatment  of  shoreline  regions  requires  special 
care:  the  water  surface  is  extrapolated  horizontally  from  the  wet  node(s)  of  a  partially 
wet  element  and  the  shoreline  defined  as  the  intersection  of  the  surface  with  the  planar 
element  topography.  This  means  that  the  shoreline  can  lie  within  an  element,  rather 
than  being  confined  to  lying  along  element  boundaries,  which  would  be  the  case  if  the 
water  surface  was  interpolated  normally  between  wet  and  dry  nodes. 

The  results  of  calibration  for  the  LISELOOD-FP  model,  using  the  kinematic  wave 
approximation  for  channel  flows  and  the  cellular  approximation  over  the  floodplain, 
are  summarised  in  table  6.1.  The  optimum  calibration  is  (by  chance)  in  the  centre  of 
the  parameter  space,  occurring  at  friction  values  we  would  expect  for  this  channel  and 
floodplain  environment.  The  model  also  shows  more  sensitivity  to  channel  friction 
than  floodplain  friction,  the  measure  of  fit  being  >80%  for  all  values  of  floodplain 
friction  when  Manning’s  n  for  the  channel  is  set  at  0.03.  The  results  of  the  calibration 
are  also  shown  in  figure  6.2  as  a  contour  plot  of  measure  of  fit  over  the  parameter 
space,  showing  the  optimum  as  a  broad  peak,  and  that  the  model  has  given  a  well 
behaved  calibration  problem.  The  best  fit  solution  is  shown  in  figure  6.3  with  the  SAR 
derived  shoreline,  showing  that  the  model  has  predicted  the  inundation  extent 
reasonably. 


nn 

0.02 

0.04 

0.06 

0.08 

0.10 

nch 

0.01 

65.2 

65.8 

65.8 

0.02 

78.5 

80.1 

80.3 

0.03 

80.9 

83.3 

83.8 

83.6 

81.9 

0.04 

79.9 

78.3 

72.1 

0.05 

79.0 

65.5 

45.1 

Table  6.1:  Calibration  for  LISFLOOD-FP  model  using  the  kinemtaic  wave 
approximation  for  channel  flow,  showing  fit  with  SAR  data  (%  correct, 
corrected  to  remove  bias)  against  floodplain  friction  (nn)  and  channel  friction 

(Hch). 


Refinements  in  the  model  process  representation  may  improve  the  results.  Firstly,  the 
diffusive  approximation  to  floodplain  flow  is  used,  and  the  results  displayed  in  figure 
6.4.  The  results  show  now  overall  improvement  in  prediction  accuracy  (a  maximum 
of  -80%),  the  model  is  slightly  more  sensitive  to  floodplain  friction  than  before,  and 
the  optimum  friction  value  has  been  shifted.  Secondly,  a  diffusive  wave 
approximation  may  be  used  in  the  channel.  Figure  6.5  shows  the  water  surface 
profiles  for  two  simulations  using  the  kinematic  and  diffusive  wave  approximations. 
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Figure  6.2:  Results  of  calibration  of  the  LISFLOOD-FP  model  using  the 
kinematic  wave  approximation.  The  black  squares  correspond  to  points  in  the 
parameter  space  for  which  simulations  were  performed,  the  surface  between 
points  is  found  from  inverse-square  distance  interpolation,  and  contoured  using 
PV-Wave  software.  The  contours  are  not  meant  to  represent  a  realistic 
interpolation,  but  are  merely  a  visualisation  tool. 
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Figure  6.3:  Best  fit  solution  from  the  LISFLOOD-FP  model  using  the  kinematic 
wave  approximation  in  the  channel  and  the  cellular  model  for  floodplain  flow. 
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Figure  6.4:  Results  of  calibration  of  the  LISFLOOD-FP  model  using  the 
kinematic  wave  approximation  in  the  channel  and  the  diffusive  model  for 

floodplain  flow. 
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Figure  6.5:  Water  surface  profiles  along  channel  for  kinematic  (solid  line)  and 
diffusive  (dotted  line)  wave  models. 
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the  kinematic  model  has  predicted  physically  unrealistic  variations  in  the  free  surface 
(induced  by  the  linkage  with  floodplain  flows),  which  are  eliminated  in  the  much 
smoother  solution  to  the  diffusive  approximation.  Nevertheless,  both  solutions  have 
the  same  overall  form,  but  with  the  diffusive  approximation  reducing  water  levels 
over  the  upstream  half  of  the  reach.  A  calibration  of  the  diffusive  channel  flow 
scheme  is  given  in  table  6.2  and  the  resulting  surface  shown  in  figure  6.6.  Figure  6.7 
shows  the  best  fit  solution  for  the  diffusive  scheme.  While  the  diffusive  scheme  has 
produced  more  realistic  predictions  of  water  depth  over  the  channel,  no  overall 
improvement  in  model  fit  is  made. 


nn 

0.04 

0.06 

0.08 

0.10 

liW 

74.2 

76.1 

77.2 

Baa 

84.2 

84.0 

0.04 

mm 

M.l 

83.7 

IIHKM 

81.3 

79.7 

Table  6.2:  Calibration  for  LISFLOOD-FP  model  using  the  diffusive  wave 
approximation,  showing  fit  with  SAR  data  (%  correct,  corrected  to  remove  bias) 
against  floodplain  friction  (nn)  and  channel  friction  (nch) 

The  results  of  the  calibration  of  the  TELEMAC-2D  model  are  given  in  table  6.3  and 
shown  in  figure  6.8.  The  model  fit  covers  a  smaller  range  than  for  LISFLOOD,  and 
the  fit  surface  has  a  more  complex  form,  although  the  higher  sensitivity  to  channel 
friction  is  still  present.  The  optimum  fit  is  similar  to  that  produced  by  the  raster  model 
using  the  kinematic  wave  approximation.  The  best  TELEMAC-2D  solution  is  shown 
in  figure  6.9  with  the  SAR  shoreline. 


nfi 
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0.04 

0.06 

0.08 

0.10 

0.12 

rich 

0.01 

75.5 

0.02 

79.0 

79.7 

0.03 

79.8 

80.2 

81.2 

81.4 

0.04 

82.6 

82.9 

83.5 

81.8 

0.05 

82.4 

82.3 

82.4 

82.6 

0.06 

82.0 

Table  6.3;  Calibration  for  TELE1VIAC-2D,  showing  fit  with  SAR  data  (% 
correct,  corrected  to  remove  bias)  against  floodplain  friction  (nn)  and  channel 

friction  (nch). 

The  raster  model  seems  to  be  capable  of  predicting  inundation  extent  reasonably  well, 
despite  the  rather  crude  assumptions  of  uniform  bed  slope,  width,  channel  depth  and 
friction.  Refinements  in  the  channel  and  floodplain  flow  representation  have  yielded 
no  improvement  in  model  results.  It  is  interesting  to  note  that  the  optimum  channel 
friction  gives  a  bankful  discharge  of  36  m^s'\  close  to  the  Environment  Agency 
estimate  of  40m^s"'.  Table  6.4  explores  the  model  sensitivity  to  channel  specification. 
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Figure  6.6:  Calibration  surface  for  the  LISFLOOD-FP  model  using  the  diffusive 
wave  approximation  in  the  channel  and  the  cellular  model  over  the  floodplain. 
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Figure  6.7:  Best  fit  solution  from  the  LISFLOOD-FP  model  using  the  diffusive 

wave  approximation. 
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Figure  6.8:  Calibration  surface  for  the  TELEMAC-2D  model. 
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Figure  6.9:  Best  fit  solution  from  the  TELEMAC-2D  model. 
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showing  that  it  is  still  possible  to  achieve  good  results  from  the  model  even  using  the 
incorrect  channel  parameters,  as  long  as  the  bankful  discharge  is  approximately 
correct.  This  may  be  indicative  of  a  certain  insensitivity  of  model  results  to  channel 
parameters,  as  long  as  the  correct  volume  of  water  is  being  routed  over  the  floodplain. 
However,  it  is  unlikely  that  this  result  will  hold  for  other  reaches  and  other  (possibly 
dynamic)  events,  and  may  be  simply  a  peculiarity  of  this  particular  model  and  data 
set. 


rich 

Width  (m) 

Depth  (m) 

Bankful  Discharge 
(m's-’) 

Fit  (%) 

0.03 

10 

3.03 

36.6 

81.3 

0.03 

20 

2.00 

36.6 

83.8 

0.03 

30 

1.57 

36.6 

83.6 

0.03 

40 

1.32 

36.6 

83.5 

0.03 

80 

0.87 

36.6 

83.1 

0.02 

20 

1.57 

36.6 

84.2 

0.04 

20 

2.00 

27.4 

78.3 

0.02 

20 

2.00 

54.9 

80.1 

Table  6.4:  LISFLOOD-FP  model  sensitivity  to  channel  specification. 


The  numerical  performance  of  the  kinematic  and  diffusive  raster  models  and  the  finite 
element  scheme  can  also  be  compared  in  terms  of  mass  balance  errors  and 
computational  efficiency,  given  in  table  6.5.  Mass  balance  errors  are  calculated  over 
each  time  step  as: 

V  -V 

Q  =Q,  _Q  — lid.  (6.1) 

Qin  is  the  imposed  upstream  flowrate,  Qom  is  the  model  downstream  flowrate,  Vt  and 
Vt-i  are  the  volumes  of  water  in  the  model  domain  at  the  current  and  previous  time 
step  and  At  is  the  model  time  step.  The  error  can  be  thought  of  in  terms  of  the  volume 
lost  or  gained  per  second  by  the  models.  While  it  is  difficult  to  come  up  with 
objective  criteria  of  what  constitute  adequate  mass  conservation  properties,  the  error 
figures  given  in  table  6.5  are  all  probably  less  than  the  error  in  the  inflow  figure  used 
to  provide  the  upstream  boundary  condition.  Given  that  the  continuity  equation  is  also 
liable  to  process  representation  errors  (it  neglects  infiltration,  runoff  from  bounding 
slopes,  rainfall  and  evaporation),  the  mass  balance  errors  found  here  probably  have  an 
insignificant  effect  on  the  predicted  inundation  extent.  The  time  taken  for  1000 
iterations,  along  with  numerical  parameters  for  the  models,  are  also  given.  The 
diffusive  raster  model  required  a  0.5s  time  step  and  10  sub  iterations  for  channel 
flows  to  achieve  stability,  compared  to  1.0s  and  no  extra  channel  sub-iterations  for  the 
kinematic  scheme.  The  small  time  step  used  is  a  result  of  instabilities  caused  mainly 
by  the  linkage  of  channel  and  floodplain  flows,  the  reduction  of  the  time  step  and  the 
use  of  sub-iterations  being  the  easiest  way  around  this  problem  without  reformulation 
the  explicit  model.  The  raster  models  were  coded  in  C-I-+  and  run  on  a  400Mhz 
Pentium  ft  processor,  and  TELEMAC-2D  in  FORTRAN  on  a  MIPS  RISC  12000 
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300MHz  processor  in  a  Silicon  Graphics  Octane  workstation.  The  figures  show  that 
the  raster  models  have  a  considerable  speed  advantage  over  the  finite  element  scheme, 
and  all  the  models  exhibit  similar  mass  balance  errors.  Given  the  size  of  the  domain 
(76  X  48  cell),  the  speed  of  the  kinematic  raster  model  is  3.5  x  10'^  s  per  cell  per  time 
step,  50  times  faster  than  the  PC-Raster  coded  model  used  in  Bates  and  De  Roo 
(2000),  with  a  speed  of  1.7  x  10""^  s  per  cell  per  time  step. 


Model 

Time 
step  (s) 

Sub-iterations 

Time  per  1000 
iterations  (s) 

Absolute 
Mass  Balance 
Error  (m^s'’) 

LISFLOOD-FP  kinematic 

1.0 

1 

12.8 

0.05 

LISFLOOD-FP  diffusive 

0.5 

10 

24.4 

0.03 

TELEMAC-2D 

2.0 

- 

169 

0.02 

Table  6.5:  Computational  performance  and  mass  balance  errors  for  the  three 

models 

As  a  final  test  of  the  performance  of  both  models,  it  is  useful  to  compare  their  results 
with  those  from  a  crude  predictor  of  flood  extent,  such  as  a  planar  water  free  surface 
height  intersected  with  the  DEM.  This  is  a  useful  safeguard  against  the  assumption 
that  the  hydraulic  models  are  performing  well,  when,  in  fact,  the  problem  of  flood 
extent  prediction  may  be  trivial.  The  planar  surface  must  be  parameterised  in  terms  of 
the  heights  at  the  upstream  and  downstream  ends,  and  several  techniques  present 
themselves.  The  heights  can  be  taken  from  water  elevations  in  the  channel,  as 
measured  at  gauging  stations,  or  from  water  elevations  taken  from  the  intersection 
between  the  SAR  shoreline  and  the  DEM.  Both  of  these  techniques  may  be  prone  to 
error  caused  by  small  local  variations  in  water  elevation  in  response  to  local  hydraulic 
conditions,  so  a  calibration  methodology  was  instead  adopted,  whereby  the  upstream 
and  downstream  elevations  of  the  planar  surface  are  adjusted  and  the  optimum  fit  with 
the  SAR  data  sought.  The  results  are  given  in  figure  6.10,  which  shows  a  definite 
maximum  (65%),  and  this  best  fit  is  shown  in  figure  6.1 1.  The  calibration  surface  for 
the  planar  predictor  shows  a  greater  sensitivity  to  downstream  height  as  the  flow  is 
less  constrained  at  the  downstream  end  due  to  the  lower  transverse  floodplain 
gradients.  Given  that  using  the  unbiased  measure  of  fit  for  classifying  the  whole  of  the 
floodplain  as  dry  is  50%,  the  results  of  the  hydraulic  models  (both  raster  and  finite 
element)  are  a  significant  improvement  over  the  planar  free  surface  predictor.  This  is 
to  be  expected  from  the  water  surface  profiles  of  figure  6.5,  which  show  a 
considerable  curvature  in  the  surface,  which  will  therefore  be  only  poorly  represented 
by  a  linear  approximation. 

6.3  Discussion 

Both  the  raster  model  and  the  finite  element  scheme  have  performed  reasonably  when 
compared  with  the  satellite  imagery,  but  it  is  unclear  whether  further  progress  can  be 
made  at  this  stage.  Errors  in  the  SAR  derived  shoreline  are  of  the  order  of  50m, 
according  to  Horritt  et  al,  2000,  who  compared  SAR  derived  flood  shorelines  against 
airphoto  data  over  two  15km  reaches  of  the  river  Thames.  If  this  error  is  reproduced 
all  along  the  shoreline,  it  is  equivalent  to  misclassifying  ~4%  of  the  domain.  This 
equates  to  an  unbiased  measure  of  fit  of  87%,  only  slightly  higher  than  the  measures 
of  fit  we  are  obtaining  with  the  inundation  models  tested  here.  This  implies  that  we 
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Figure  6.10:  Calibration  of  planar  free  surface  approximation,  showing  measure 
of  fit  as  a  function  of  upstream  and  downstream  water  elevations. 
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Figure  6.11:  Best  fit  solution  using  the  planar  free  surface  approximation. 
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are  achieving  a  similar  measure  of  fit  between  model  and  SAR  data  as  between  SAR 
data  and  the  real  flood  shoreline,  and  any  further  improvement  in  model  performance 
is  not  possible  with  this  data  set.  Both  the  raster  and  the  finite  element  models  achieve 
a  similar  measure  of  fit,  but  to  distinguish  further  between  them  requires  more  precise 
validation  data.  It  must  be  remembered  that  the  optimum  model  formulation  and 
calibration  is  not  only  dependent  on  the  validation  data,  but  the  objective  function 
used  in  the  comparison  with  model  predictions.  The  results  of  inappropriate  process 
representation  have  thus  been  masked  by  uncertainty  in  the  validation  data,  and  at  this 
level  of  uncertainty,  a  relatively  crude  process  representation  is  enough  to  reproduce 
the  observed  inundation  extent. 

It  should  be  stressed  that  this  research  has  focussed  on  the  use  of  hydraulic  models  as 
inundation  predictors,  and  they  have  been  tested  as  such  against  the  SAR  derived 
shoreline  data.  There  will  be  only  a  weak  link  between  inundation  extent  and 
floodplain  and  channel  hydraulics,  so  this  research  only  partially  validates  the 
hydraulic  representations  used  in  the  model.  This  is  demonstrated  particularly  well  by 
the  raster  model’s  robustness  with  respect  to  channel  specification.  Very  different 
channel  widths  and  depths,  with  the  constraint  that  they  should  approximately 
reproduce  bankful  discharge,  can  produce  similar  inundation  patterns  and  measures  of 
fit  when  compared  to  the  SAR  data,  but  with  very  different  hydraulic  conditions 
operating  in  the  channel.  This  property  is  also  demonstrated  by  the  similarity  of  the 
raster  models’  predictions  of  flood  extent  using  the  kinematic  and  diffusive  wave 
approximation.  The  unphysical  variations  in  the  free  surface  over  the  channel 
predicted  by  the  kinematic  model  are  not  reflected  in  the  flood  shoreline,  perhaps 
because  these  variations  are  in  some  sense  damped  out  in  the  far  flow  field  near  the 
shoreline.  This  is  encouraging  in  terms  of  inundation  prediction,  as  it  appears  that  in 
this  case,  as  long  as  bankful  discharge  is  approximately  correct,  the  inundated  area 
will  also  be  approximately  correct.  Looking  at  the  inverse  problem,  however,  we  see 
that  inundation  extent  may  give  us  very  little  information  about  channel  flows. 

This  study  has  been  limited  to  a  single  test  site,  and  further  applications  to  different 
sites  are  required  to  verify  or  disprove  the  results  obtained  here.  There  are  a  number 
of  factors  that  may  affect  prediction  accuracy  when  these  modelling  strategies  are 
extended  to  other  reaches.  Inundation  extent  is  very  sensitive  to  topography,  so  the 
details  of  floodplain  topography  will  affect  the  sensitivity  to  model  formulation  and 
calibration.  Dynamic  effects  may  also  become  important,  especially  over  longer 
reaches.  The  Idnematic  wave  velocity  for  the  channel  used  here  will  be  approximately 
1.5  ms"’  (Chow,  1988,  p284),  and  the  surface  wave  celerity  approximately  4.5  ms’’,  so 
kinematic  effects  should  propagate  along  the  4km  reach  in  '-1  hour.  The  maximum 
rate  of  change  of  the  hydrograph  for  this  event  is  1.3  m^s‘'hT*,  so  dynamic  effects  are 
unlikely  to  be  important  for  this  reach.  Longer  reaches  will  respond  more  slowly  and 
with  more  complexity  to  the  dynamic  behaviour  of  the  input  hydrograph,  and  this  may 
become  a  useful  diagnostic  tool  for  assessing  model  performance.  For  example,  we 
have  seen  how  the  raster  based  model  is  fairly  robust  with  respect  to  channel 
specification,  as  long  as  bankful  discharge  is  reproduced  reasonably  well.  Varying 
channel  depth  does,  however,  also  affect  kinematic  and  diffusive  wave  velocities,  and 
so  changing  the  channel  properties  may  have  much  more  of  an  effect  for  dynamic 
simulations  than  the  steady  state  solutions  developed  here. 
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The  application  of  hydraulic  models  to  inundation  prediction  is  still  reliant  on  the 
calibration  process,  which  tends  to  obscure  model  validation  issues.  The  calibration 
process  is  model  dependent:  the  TELEMAC-2D  model  produces  the  most  complex 
calibration  surface,  while  that  for  the  LISFLOOD-FP  model  with  kinematic 
approximation  for  channel  flows  is  relatively  simple,  with  the  diffusive  model 
somewhere  in  between.  With  all  the  models  giving  roughly  the  same  level  of  fit  at  the 
optimum  calibration  (the  optimum  being  model  dependent),  this  is  a  point  in  favour  of 
the  simpler  model.  Given  that  currently  we  can  only  achieve  closure  in  terms  of 
friction  parameterisation  through  a  calibration  procedure,  the  simpler  calibration 
properties  of  the  raster  model  are  a  definite  advantage. 

6.4  Conclusions 

All  the  models  tested  here  (raster  kinematic,  raster  diffusive  and  finite  element)  have 
performed  to  a  similar  level  of  accuracy,  classifying  approximately  84%  of  the  model 
domain  correctly  when  compared  to  SAR  derived  shoreline  data.  However,  the 
validation  data  used  here  is  insufficiently  accurate  to  distinguish  between  the  model 
formulations,  and  the  issue  of  model  validation  is  further  clouded  by  the  calibration 
process  necessary  due  to  the  lack  of  friction  parameterisation  data.  The  remotely 
sensed  data  has  proved  invaluable,  however,  in  the  calibration  process,  and  has 
reduced  the  equifinality  problem  inherent  in  calibrating  distributed  models  with  point 
hydrometric  data.  Given  the  likely  accuracy  of  the  validation  data,  and  the  increasing 
complexity  of  the  calibration  process  for  models  representing  more  complex 
processes,  the  simple  raster  based  model  using  a  kinematic  wave  approximation  over 
the  channel  is  the  simplest  (and  fastest)  model  to  use  and  adequate  for  inundation 
prediction  over  this  reach. 


7.  The  impact  of  varying  levels  of  bathymetric  data  and  mesh 
resolutions  on  the  LISFLOOD-FP  model  predictions. 

Section  6  established  that  even  for  best  available  data  sets  it  may  not  be  possible  to 
discriminate  between  simple  and  complex  hydraulic  models.  Given  the  relatively 
good  performance  of  the  LISFLOOD  model,  its  computational  efficiency  and  raster 
basis  enable  a  more  rigorous  examination  of  data  and  resolution  requirements  than 
would  have  be  possible  with  the  TELEMAC-2D  model.  This  work  therefore  has  the 
potential  to  clarify  a  number  of  areas  concerning  topographic  data  provision  for 
hydraulic  models  that  can  then  be  translated  to  more  complex  codes. 

To  achieve  this,  models  of  resolution  varying  from  10m  to  1000m  were  built  for  a 
60km  reach  of  the  River  Severn,  UK  and  predictions  compared  with  satellite  SAR 
observations  of  inundated  area  and  ground  measurements  of  floodwave  travel  times, 
with  a  calibration  strategy  being  used  to  determine  channel  friction  coefficients.  The 
optimum  calibration  was  found  to  be  stable  with  respect  to  changes  in  scale  when  the 
model  was  calibrated  against  the  observed  inundated  area,  the  model  reaching 
maximum  performance  at  a  resolution  of  100m,  after  which  no  improvement  was  seen 
with  increasing  resolution.  Projecting  predicted  water  levels  onto  a  high  resolution 
DEM  improved  performance  further,  and  a  resolution  of  500m  proved  adequate  for 
predicting  water  levels.  Predicted  floodwave  travel  times  were,  however,  strongly 
dependent  on  model  resolution,  and  water  storage  in  low  lying  floodplain  areas  near 
the  channel  was  identified  as  an  important  mechanism  affecting  wave  propagation 
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velocity.  A  near  channel  floodplain  storage  version  of  the  LISFLOOD-FP  model  was 
shown  to  be  much  more  stable  with  respect  to  changes  in  scale  when  the  model  was 
calibrated  against  floodwave  travel  times,  and  shown  to  represent  the  retardation  of 
the  floodwave  caused  by  water  storage  near  the  channel.  The  LISFLOOD-FP  model 
could  not,  in  this  instance,  be  calibrated  to  simultaneously  give  both  acceptable  travel 
times  and  inundated  areas. 


8.  Future  R&D  opportunities 

This  project  has  sought  to  develop  a  suite  of  computational  hydraulic  models  for  high- 
resolution  flow  prediction  at  the  reach  scale  (10-60km)  that  directly  addresses 
potential  COE  Research  and  Development  agendas.  Specifically,  we  have: 

•  Developed  a  suite  of  models  of  varying  complexity  for  long  reach,  high-resolution 
river  flow  prediction. 

•  Developed  through  GIS  technologies  the  integration  to  remote  sensing  data  sets 
capable  of  parameterizing  such  models  and  examined  data  assimilation, 
redundancy  and  scaling  issues. 

•  Developed  novel  means  of  validating  hydraulic  models  using  newly  available  data 
sets  from  satellite  and  airborne  sensing  platforms. 

Future  opportunity  :  With  the  further  advance  of  remote  sensing  technologies  and 
high  performance  computing  considerable  potential  now  exists  for  computational 
hydraulic  modelling  at  all  scales  up  to  and  including  the  basin  scale  (lOO’s  of  km). 
This  has  the  potential  to  allow  extension  of  the  modelling  techniques  described  in  this 
report  into  the  areas  of: 

Forecasting 

•  Real  time  forecasting 

•  Linkages  to  remote  sensing  driven  snow-melt  forecasting  models 

•  Hydraulic  impacts  of  climate  and  land  use  change 

Design 

•  Soft  engineering  design  in  respect  of  land  use,  habitat  specifications 

•  Development  of  maintenance  schedules  and  impact  assessments. 

Management 

•  Linking  hydrologic  and  hydraulic  models  for  Integrated  Basin  Management 

•  Floodplain  management  and  planning,  including  integration  of  model  outputs  with 
socio-economic  data  sets  to  identify  at  risk  populations 

•  Wetland  restoration  and  management 
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